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The  principal  objective  of  this  research  is  to  extend  the  results 
of  R.  E.  Kalman's  original  analysis  of  the  inverse  optimal  control  prob- 
lem in  order  to  create  a  design  tool  which  allows  the  designer  to  employ 
the  powerful  and  sophisticated  techniques  of  optimal  control  theory 
when  the  correct  choice  of  a  performance  index  is  not  apparent.   Specif- 
ically, methods  are  developed  for  the  design  of  optimal  linear  regulators 
(both  continuous  and  discrete)  which  satisfy  classical  performance 
specifications  in  addition  to  the  minimality  of  a  quadratic  functional. 
This  permits  the  great  computing  power  of  optimal  control  theory  to  be 
brought  to  bear  on  problems  which  previously  could  be  accommodated  only 
by  the  cut-and-try  methods  of  classical  design  schemes.   In  passing, 
insights  into  the  basic  nature  of  optimal  regulators  are  developed  in 
terms  of  classical  concepts,  unifying  some  important  notions  in  clas- 
sical and  modern  control  theory. 

The  contributions  of  this  research  can  be  summarized  as  follows: 
1.   A  complete  theoretical  investigation  of  the  inverse 
optimal  linear  regulator  problem  for  quadratic  per- 
formance indices  with  positive  semidefinite  state 


weighting  matrices  and  scalar  input  systems  is  reviewed 
and  preliminary  results  for  the  case  where  the  weighting 
matrix  is  allowed  to  be  sign  indefinite  are  given. 

2.  The  equivalence  of  performance  indices  for  scalar  input 
linear  quadratic  loss  problems  is  resolved  and  a  proce- 
dure for  generating  the  entire  equivalence  class  of  cost 
functions  equivalent  to  a  given  performance  index  is 
developed. 

3.  Practical  numerical  methods  for  determination  of  system 
optimality  and  computation  of  solutions  to  the  scalar 
input  inverse  problem  are  discussed. 

4.  Some  of  the  effects  of  specific  elements  of  the  perform- 
ance index  on  optimal  system  performance  and  pole  loca- 
tions are  determined. 

5.  The  problem  of  designing  optimal  systems  to  meet  classical 
performance  specifications  is  encountered  and  some  defin- 
itive results  obtained. 

6.  A  solution  to  the  problem  of  designing  a  sampled-data 
controller  to  approximate  the  performance  of  continuous 
controller  is  given. 


CHAPTER  I 


INTRODUCTION 


1. 1   Background 

In  1964  R.  E.  Kalman  published  a  paper  [Kl]   which  has  come  to 
have  considerable  impact  on  the  theory  of  optimal  control.   It  dealt 
not  with  the  conventional  optimal  control  problem  of  computing  a  system 
trajectory  which  extremizes  a  specified  performance  index,  but  rather 
with  the  "inverse"  problem  of  determining  what  performance  indices,  if 
any,  are  extremized  by  a  specified  control.   Kalman  restricted  his 
attention  to  the  case  where  the  system  is  linear  and  the  integral  per- 
formance index  is  quadratic  in  the  states  and  control.   No  immediate 
direct  application  was  made  of  Kalman' s  results,  but  they  were  of  great 
value  in  the  application  and  interpretation  of  his  earlier  analysis 
[K2]  of  what  is  now  called  the  "linear  quadratic  loss"  optimization 
problem. 

Quadratic  functionals  have  long  been  studied  in  the  calculus 
of  variations  [Gl]  but  not  until  Kalman  was  the  minimization  of  a  func- 
tional quadratic  in  the  states  and  input  of  a  linear  dynamical  system 
considered  in  the  context  of  an  optimal  control  problem.   Other  authors, 
notably  S.  S.  L.  Chang  TCI]  and  Newton,  Gould,  and  Kaiser  [Nl] , 
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betic group  of  the  specific  reference. 


had  since  the  middle  1950'  s  been  occupied  with  the  "analytical  design"  of 
linear  control  systems  through  the  use  of  a  performance  index  quadratic 
in  the  error  (i.e.,  the  difference  between  actual  and  desired  system 
responses  to  a  specified  input).   Employment  of  Parseval's  Theorem  and 
spectral  factorization  [Nl]  resulted  in  a  solution  strongly  related  to 
the  optimal  linear  filter  of  Wiener  [Wl]  of  almost  fifteen  years  earlier. 
The  solutions  for  the  ISE  (integral  of  the  square  of  the  error)  problems 
were  very  cumbersome  to  compute  and  seldom  applied  to  systems  of 
greater  than  third  order. 

In  the  early  1960's  the  fledgling  field  of  optimal  control 
theory  underwent  a  metamorphosis.   Difficult  aerospace  problems  had 
arisen  which  could  only  be  conveniently  approached  as  optimal  control 
problems.   Rapidly  there  began  appearing  almost  as  many  different 
optimal  control  formulations  as  there  were  proponents.   Superficially, 
optimal  control  theory  appeared  to  be  a  panacea;  however,  formidable 
problems  in  computing  and  often  even  greater  difficulties  in  the 
implementation  of  some  optimal  designs  limited  their  utility.   In 
general,  the  computation  of  optimal  controls  requires  the  use  of 
iterative  algorithms  which  converge,  at  best,  slowly.   Further,  optimal 
controls  are  usually  open-loop  in  nature. 

The  linear  quadratic  loss  problem  suffers  from  neither  of 
these  difficulties:  the  solution  is  numerically  straightforward  and 
always  results  in  a  linear  feedback  control  law  which  is  easily 
implemented;  in  addition  both  continuous  and  sampled-data  controllers 
can  be  accommodated.   If  the  system  is  to  operate  in  a  noisy  environment, 
an  extension  of  the  quadratic  loss  formulation,  the  Kalman-Bucy  filter, 
can  easily  be  included  in  the  design. 


Unfortunately,  performance  indices  as  a  whole  seldom  relate 
well  to  system  design  requirements  which  are  given  as  classical  time 
and  frequency  .domain  specifications.   In  this  regard,  the  quadratic 
loss  problem  suffers  as  well. 

The  principal  objective  of  this  research  is  to  extend  those 
original  results  of  Kalman  in  order  to  create  a  design  tool  which 
allows  the  designer  to  employ  the  powerful  and  sophisticated  techniques 
of  optimal  control  theory  when  the  correct  choice  of  a  performance 
index  is  not  apparent.   Specifically,  methods  are  developed  for  the 
design  of  optimal  linear  regulators  (both  continuous  and  discrete) 
which  satisfy  classical  performance  specifications  in  addition  to  the 
minimality  of  a  quadratic  functional.   This  would  permit  the  great 
computing  power  of  optimal  control  theory  to  be  brought  to  bear  on 
problems  which  previously  could  be  accommodated  only  by  the  cut-and-try 
methods  of  classical  design  schemes.   In  passing,  insights  into  the 
basic  nature  of  optimal  regulators  are  developed  in  terms  of  classical 
concepts,  unifying  some  important  notions  in  classical  and  modern 
control  theory. 

The  original  impetus  for  this  investigation  was  the  requirement 
to  develop  a  method  for  designing  a  digital  controller  to  replace  an 
existing  continuous  compensator  without  significantly  affecting  system 
performance.   It  was  believed  that  determination  of  a  performance  index 
minimized  by  the  continuous  system  would  allow  for  performance  invariant 
design  by  computing  the  optimal  control  law  for  the  sampled-data  ver- 
sions of  the  continuous  system  and  performance  index.   In  Chapter  IV 
this  scheme  is  discussed  as  a  solution  to  a  problem  which  will 


undoubtedly  occur  with  increasing  frequency  as  digital  controllers 
become  more  commonplace.   It  was  in  the  course  of  this  investigation 
that  the  potential  of  the  inverse  problem  for  regulator  design  of  wider 
latitude  was  discovered. 

1.2   Survey  of  Previous  Work 

Perhaps  the  first  attempts  at  relating  optimal  control  theory 
and  classical  design  were  early  (c.  1950)  investigations  of  so-called 
"standard  forms"  (e.g.,  [G2]).   The  object  was  to  tabulate  forms  of 
closed-loop  system  transfer  functions  which  were  optimal  with  respect 
to  a  specified  performance  measure  (for  instance,  ISE)  and  input. 
This  approach  was  rather  restrictive  and  difficult  to  apply  to  problems 
of  interest.   Its  impact  was  nonetheless  considerable  and  a  recent 
paper  by  Rugh  [Rl]  indicates  that  the  linear  quadratic  loss  problem 
has  much  in  common  with  these  early  results. 

Shortly  after  the  linear  quadratic  loss  problem  became  widely 
known,  the  task  of  selecting  performance  indices  which  result  in 
optimal  systems  with  desired  characteristics  was  attacked  primarily  on 
an  experimental  basis.   The  hope  was  that  massive  experience  and  some 
insight  into  the  mechanics  of  computing  the  optimal  control  laws  would 
lead  to  guidelines  for  the  choice  of  a  quadratic  performance  index. 
One  such  procedure,  developed  by  Tyler  and  Tuteur  [Tl]  consisted  of 
computing  root-locus  plots  as  functions  of  weighting  terms  in  the 
performance  index  and  choosing  a  suitable  compromise. 

Another  study  of  optimal  systems  in  classical  terms  was  made 
in  1965  by  R.  J.  Leake  [LI].   Leake,  in  developing  a  computational 
scheme  for  solutions  to  the  linear  quadratic  loss  problem,  shows  how 


an  estimate  of  optimal  system  bandwidth  may  be  made  without  computing 
the  optimal  control  and  comments  on  an  intriguing  geometric  interpreta- 
tion (in  the  complex  plane)  of  Kalman' s  criteria  for  optimality. 

An  important  step  in  the  analysis  that  follows  will  be  the 
determination  of  when  two  different  quadratic  performance  indices 
applied  to  the  same  plant  will  result  in  identical  optimal  control  laws. 
Kalman' s  paper  on  the  inverse  problem,  by  its  nature,  encounters  the 
problem  theoretically  but  does  not  consider  it  in  detail.   A  partial 
practical  solution  to  the  equivalence  problem  appeared  simultaneously 
with  Kalman' s  paper.   Wonham  and  Johnson  observed  in  their  paper  on 
the  linear  quadratic  problem  with  a  bounded  control  [W2]  that  an  arbi- 
trary weighting  matrix  in  the  performance  index  could  be  replaced  by 
a  diagonal  one  which  will  result  in  the  same  optimal  control.   Their 
result  arises  from  the  observation  that  when  the  system  is  expressed 
in  a  canonical  form  (companion  matrix  form  [R2])  the  performance  index 
may  be  reduced  by  repeated  integration  by  parts  to  a  diagonal  form; 
they  fail  to  realize,  however,  that  the  diagonalized  version  may  no 
longer  possess  a  solution.   Kalman  and  Englar  later  recognize  [K3, 
p.  304-306]  that,  in  the  same  canonical  form  as  used  by  Wonham  and 
Johnson,  certain  terms  may  be  discarded  and  equivalence  maintained. 

Beyond  these  early  results  and  an  occasional  rediscovery  of 
them  (e.g.,  Kreindler  and  Hedrick  [K4])  the  study  of  equivalent  quad- 
ratic performance  indices  has  remained  dormant. 


In  summary: 

1.  A  great  deal  of  experimental  work  has  been  done  in  hopes  of 
relating  optimal  design  to  classical  criteria  and  providing 
intuition  into  a  procedure  for  choosing  performance  indices. 
Very  little  basic  theoretical  research  has  been  invested 

in  this  area. 

2.  The  practical  problems  of  determining  the  optimality  of  an 
actual  system  and  the  computation  of  a  pei'formance  index 
have  not  been  considered. 

3.  Some  elementary  equivalence  relations  have  been  developed, 
but  the  difficulties  arising  from  naively  applying  them 
have  been  generally  overlooked.   Nor  have  the  advantages  of 
one  equivalent  form  over  another  been  explored. 

4.  No  real  effort  has  been  made  to  use  the  study  of  the  inverse 
problem  as  a  foundation  for  linear  regulator  design. 

This  work  begins  with  a  review  of  the  linear  plant,  quadratic  cost 
variational  problem.  The  conditions  for  the  solution  to  exist  are  dis- 
cussed along  with  techniques  for  computing  the  optimal  control. 

The  inverse  problem  for  a  linear  plant  and  quadratic  performance 
index  is  next  considered.   The  conditions  for  optimality  are  developed, 
the'  meaning  of  optimal  control  is  studied  in  terms  of  classical  criteria 
and  the  equivalence  of  performance  indices  is  resolved. 

The  fourth  chapter  develops  computational  procedures  for  the 
determination  of  solutions  to  the  inverse  problem  and  goes  on  to  con- 
sider practical  techniques  for  the  design  of  optimal  regulators  which 
meet  classical  performance  specifications. 


CHAPTER  II 
THE  LINEAR,  QUADRATIC  COST  OPTIMIZATION  PROBLEM 

2. 1   Introduction 

As  a  first  step  toward  the  establishment  of  viable  design 
techniques  based  on  the  linear,  quadratic  cost  optimal  control  problem, 
the  optimization  problem  itself  must  be  reviewed  in  some  detail.   The 
great  power,  latitude  and  computational  elegance  of  this  optimal  control 
formulation  account  in  large  part  for  its  popularity  as  an  object  of 
study.   In  the  investigation  that  follows  its  limitations  will  also 
have  considerable  impact. 

This  chapter  will  first  define  what  is  meant  by  the  linear 
plant,  quadratic  cost  optimization  problem.   A  method  for  the  compu- 
tation of  solutions  which  also  provides  a  necessary  and  sufficient 
test  for  existence  will  then  be  considered.   The  chapter  closes  with 
a  specific  review  of  the  particular  subproblem  which  will  be  of 
principal  interest  in  the  remainder  of  this  work,  the  time-invariant 
optimal  linear  regulator. 


2.2   The  Linear  System,  Quadratic  Cost  nroblem 

Consider  an  n-dimensional  linear  system  with  state  feedback 

described  by  the  state  equation 

dx 

^=Fx+Gu  (1) 

T 
where  u  =  -K  (t)x 

and  a  companion  performance  index 

t 


1  /»    T      T  T  T 

J  =  -  J   (xQx+u  Ru)dt,   Q  =  Q   and  R  =  R  (2) 


where  u  is  m-dimensional ;  F,  G,  Q,  and  R  are  real,  constant  matrices; 
and  K  is  a  real  (possibly  time-varying)  matrix;  all  are  of  appropriate 
dimensions. 

It  will  be  shown  later  (Appendix  A)  that  several  other 
performance  indices  can  be  accommodated  in  the  framework  established 
for  this  particular  one.   The  —  preceding  the  integral  in  (2)  is  for 
algebraic  simplicity  in  the  present  discussion  and  has  no  effect  on 
the  result  of  the  optimization  process;  it  will  occasionally  be 
deleted  in  the  sequel  without  comment. 

Taken  together  there  are  a  variety  of  problems  inferred  by 
(1)  and  (2).   The  most  obvious  is  the  optimal  control  problem  [K2] . 
i.e. ,  compute  the  control  law  K(t),  if  any,  which  minimizes  the 
performance  index  (2).   A  second,  somewhat  more  obscure  problem,  is 
the  so-called  "inverse  problem"  of  optimal  control  theory  [K1.E1]. 
The  "inverse  problem"  seeks  to  determine  what,  if  any,  performance 
index  is  minimized  by  a  specific  feedback  control  law.   In  the  fol- 
lowing, both  of  these  topics  will  play  an  important  role  in  uncover- 
ing the  nature  of  optimal  systems. 


2.3   The  Optimal  Control  Problem 


The  Euler-Lagrange  Equations 


The  necessary  conditions  for  a  control  u  to  minimize  the 
performance  index  (2)  is  that  the  system  equations  (1)  and  the 
celebrated  Euler-Lagrange  equations  [Gl]  be  simultaneously  satisfied, 
i.e., 


•      dH 

x  =  "  5x« 


X(tf)  =  0 


dH 

dTT 


(3) 
(4) 


where  H  is  the  Hamiltonian 


1   T_    IT, 


H=-xQx  +  -uRu  +  \    (Fx  +  Gu) 


If  R  is  non-singular,  equation  (4)  (sometimes  called  the  "stationarity 
condition"  [B3])   determines  a  candidate  for  the  optimal  control,  i.e. , 


u(t)  =  -R_1GTX(t) 


(5) 


The  simultaneous  solution  of  equations  (1)  and  (3)  with  the  substi- 
tution of  (5)  results  in  the  two-point  boundary  value  problem: 


_d_ 

dt 


-lT-ti 1 

-GR   G  n 


x(t  )  =  x   (given) , 


\(tf)  =  0 


(6) 
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The  Riccati  Equation 

The  optimal  control  problem  for  a  linear  system  subject  to 
a  quadratic  performance  index   has  been  studied  in  great  detail  over 
the  past  several  years  (e.g.,  [K2] ,  [K3])  and  several  eloquent  solu- 
tions for  the  necessary  conditions  have  been  offered  [Bl].   Most  of 
these  solutions  differ  only  in  the  approach  taken  to  solve  the  two- 
point  boundary  value  problem  (6);  for  the  purpose  of  this  chapter  the 
so-called  "sweep  method"  [Gl]  is  most  illustrative  and  will  be  pre- 
sented by  way  of  review. 

Hypothesize  a  matrix  P(t,t  )  such  that 

X(t)  =  P(t,tf)x(t),  (7) 

then  P(t,t  )  would  in  effect  provide  a  boundary  value  which  is 
"swept  back"  in  time  to  the  initial  time  and  the  initial  value  for 
X  would  simply  be 

X(tQ)  =  P(to,tf)xo, 

hence,  the  separation  of  the  boundary  values  would  be  resolved. 

By  requiring  that  the  Euler-Lagrange  equations  be  satisfied,  it  may 

be  shown  that 

/HO  T  —1  T        N 

!  -r-  +  PF  +  F  P  -  PGR   G  P  +  Q  x  =  0 . 
\dt  J 

Since  x(t)  is  arbitrary,  P(t,t  )  must  satisfy 

d£  =  _pF  _  FTp  +  PGr"1gTP  -  Q,  P(t  ,t  )  =  0  (8) 

dt  I   J 

which  is  a  matrix  version  of  the  familiar  Riccati  equation  [Dl] . 


1This  problem  is  often  referred  to  as  the  "optimal  linear 
regulator"  problem  in  the  literature;  however,  this  designation  will 
be  reserved  for  the  infinite  final  time  case  here. 
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At  first  it  may  seem  rather  surprising  that  the  solution  to 
a  2n  order  two-point  boundary  value  problem  (6)  can  be  obtained  from 
the  solution  of  an  nth  order  non-linear  equation  (8),  but  recalling  that 
the  scalar  Riccati  equation  was  solved  in  elementary  differential 
equation  theory  with  the  aid  of  a  2nd  order  linear  differential  equation 
relates  it  to  a  familiar  problem. 

Finally,  substitution  of  (7)  into  (5)  leads  to  a  candidate  for 
the  optimal  control  (a  control  which  satisfies  the  necessary  condi- 
tions) in  the  form  of  a  feedback  control  law,  i.e. , 

u  =  -KT(t)x,  where  K(t)  =  P(t,tf)GR_1  (9) 

and  P  is  the  solution  to  the  Riccati  equation  (8). 

Before  proceeding  any  further,  another  necessary  condition  is 
immediately  available  which  has  been  previously  ignored  in  this  analysis. 
From  observation  of  the  performance  index  (2) ,  it  is  obvious  that  R 
must  be  positive  semidef inite;  otherwise,  a  control  could  be  hypoth- 
esized with  sufficient  high  frequency  content  to  have  negligible  effect 
on  the  states  of  the  system,  hence  allowing  the  integral  (2)  to  become 
unbounded  (negatively).   The  requirement  in  the  Euler-Lagrange  equa- 
tions that  R  be  non-singular  further  constrains  R  to  be  positive  definite. 
This  assumption  is  equivalent  to  the  "strengthened  Legendre-Clebsche" 
necessary   condition  of  the  calculus  of  variations;  therefore,  the 
requirement  that  R  be  positive  definite  is  sometimes  referred  to  by 
that  nomenclature. 

In  summary,  a  solution  to  the  matrix  Riccati  equation  (8) 
specifies,  in  the  form  of  equations  (9),  a  control  which  satisfies 
the  necessary  conditions  (Euler-Lagrange  equations)  for  optimality; 
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nil  that  remains  is  to  investigate  sufficient  conditions  for  the 
optimal ity  of  (9).   The  existence  of  a  solution  to  the  Riccati  equation 
and  an  additional  sufficient  condition  are  related  to  the  non-existence 
of  conjugate  points  which  are  defined  below. 


The  Conjugate  Point  Condition  [Gl] 

Definition  2. 1 

A  point  t.  is  called  a  conjugate  point  to 
tf  if  there  exists  a  solution  to  the  Euler-Lagrange 
equations  (6)  with  the  boundary  conditions 

x(tj  =  x(t  )  =  0,  where  t  £  t  <  t  , 
1       f  o    1    f 

which  is  not  zero  everyv/here  on  the  interval. 
Graphically,  a  conjugate  point  can  be  illustrated  as  in  Figure  2.1 
(cf.  Figure  1  in  Reference  [B3]). 


x(t) 


Figure  2.1  Conjugate  Trajectory 


Lemma  2. 1  [B3] 

When  a  conjugate  point  exists,  both  the 
trivial  solution  (Path  1  in  Figure  2.1)  and  the 
conjugate  trajectory  (Path  2  in  Figure  2.1) 
result  in  a  zero  value  for  the  integral 

n     T       T 

(x  Qx  +  u  Ru)dt. 
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Proof 

Consider  the  identically  zero  integral  on  the  conjugate  path 


*f 


T  .  „-l„T, 


J    X  (x  -  Fx  +  GR  "G'X)dt  =  0 

J*   \Txdt  -  J*   XT(Fx  -  GR_1GTX)dt  =  0. 


Applying  integration  by  parts  to  the  first  term  results  in 


*f    *f 


T     ,T,      „r,-l«T 


XTx  L   -  JJ.   XTx  +  X  (Fx  - 


GR   G  X)dt  =  0. 


The  first  terra  is  clearly  zero  and  with  the  substitution  for  X  from 
the  Euler-Lagrange  equations  (6)  and  for  u  from  (5) 

*« 

P    (xTQx  +  uTRu)dt  =  0, 
*1 

which  completes  the  proof. 

When  a  conjugate  point  does  exist  at  t  ,  X(t  )  must  be  non-zero 
to  be  distinct  from  the  identically  zero  solution  to  (6).   If  the  system 
is  controllable,  it  must  respond  to  the  input  (5)  resulting  from 
X(t  )  4   0;  hence  by  continuity  the  conjugate  path  (X-and  x)  must'  be  non- 
zero for  some  finite  interval  within  (t   t  ).   Then  the  conjugate  path 
cannot  be  a  duplicate  of  the  identically  zero  path  with  the  addition 
of  a  (non-zero)  discontinuity;  thus  the  presence  of  a  conjugate  point 
represents  more  than  merely  misbehavior  at  a  single  point.   The  occur- 
rence of  a  conjugate  point  coincident  with  t'  is  obv.ioiislyt'disconcerting 
because  it  indicates  that  there  are  at  least  two  candidates. (that' 
satisfy  the  necessary  conditions)  for  an  optimal  control  which  lead 
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to  entirely  different  trajectories  at  the  same  (zero)  cost.   When  a 
conjugate  point  occurs  at  t  ?!  t  ,  the  results  are  equally  catastrophic 
but  this  case  does  not  lend  itself  as  well  to  heuristic  interpretation. 

Theorem  2.  1 

For  a  linear,  completely  controllable  system  subject 

to  a  quadratic  performance  index  (2)  with  R  positive 

definite,  a  control  u*  satisfying  the  Euler-Lagrange 

equations  is  globally  optimal  if  and  only  if  there  exist 

no  conjugate  points  on  the  interval  (t  ,t_). 

o   I 

A  discussion  of  this  theorem  and  its  proof  in  a  somewhat  dif- 
ferent context  can  be  found  in  Breakwell  and  Ho  [B3] .   The  requirement 
that  conjugate  points  be  non-existent  is  generally  referred  to  as  the 
"conjugate  point  condition"  or  occasionally  in  the  classical  calculus 
of  variations  as  the  Jacobi  condition  [Gl]. 

Complete  controllability  requires  that  there  exists  an  input, 
u(t) ,  which  will  drive  the  system 

x  =  Fx  +  Gu(t)  ,     x(t  )  =  x 
o     o 

from  any  initial  state,  x  ,  to  the  origin  within  an  arbitrary  time 
interval;  this  can  be  shown  to  be  equivalent  to  requiring  that  the 
matrix  [G,FG, ... ,F  ~  G]  have  full  rank  [K5].   Controllability  is  not 
actually  a  severe  criterion  in  a  practical  sense.   A  plant  which  is 
not  completely  controllable  can  be  transformed  into  two  canonical 
subsystems,  one  containing  the  completely  controllable  part  and  the 
other  subsystem  containing  the  remainder  of  the  plant  dynamics  [K5] . 
Hence,  a  given  design  problem  can  be  thought  of  as  two  related  designs 
on  the  canonical  subsystems  and  it  is  only  necessary  to  be  certain  that 
the  non-controllable  part  is  stable  and/or  not  reflected  in  the 
performance  index. 
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Now  the  solution  to  the  linear,  quadratic  cost  optimization 
problem  is  essentially  complete;  it  is  merely  necessary  to  solve  the 
Riccati  equation  for  a  candidate  control  law  and  to  insure  that  the 
conjugate  point  condition  is  satisfied.   Although  there  appears  to  be 
no  convenient  way  to  test  for  the  presence  of  conjugate  points,  the 
following  theorem  demonstrates  that  a  test  for  the  conjugate  point 
condition  is  actually  implicit  in  the  solution  of  the  Riccati  equation. 

Theorem  2. 2 

The  solution  to  the  Riccati  equation  (8)  for  the 
optimization  problem  specified  by  (1)  (completely  con- 
trollable) and  (2)  (with  R  positive  definite)  fails  to 

exist  (becomes  unbounded)  at  t,  ,  t  st,  <t_,  if  and  only 
....  .  ,     1.   o    1   I 

if  t   is  a  conjugate  point  to  t  . 

A  formal  proof  can  be  found  in  Lee  [L2]  which  is  very  much  in 
the  spirit  of  the  following  heuristic  justification. 

From  the  definition  of  a  conjugate  point,  x(t  )  =  0;  however, 
\(t  )  ^  0,  since  this  could  only  come  about  in  the  trivial  solution. 
Recall  that  the  solution  to  the  Riccati  equation  was  defined  by  (7)  as 

X(t)  =  P(t,tf)x(t). 

Then  as  t  approaches  t  ,  x  approaches  zero  but  \   remains  non-zero; 
hence,  ||P(t)||  must  correspondingly  increase  and  finally  become 
unbounded  at  t  =  t  . 

Thus  far  the  conjugate  point  condition  has  been  presented  in 
a  rather  esoteric  format  with  little  physical  interpretation.   Again 
consideration  of  the  general  case  is  difficult  but  the  occurrence  of 
a  conjugate  point  at  the  initial  time  (t  )  has  an  interesting 
interpretation.   The  absence  of  a  conjugate  point  at  t   (in  the  linear 
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case)  insures  that  no  two  system  trajectories  which  satisfy  the  Euler- 

2 

Lagrange  equations  (extremals)  will  ever  intersect,   that  is,  there  will 

never  exist  two  distinct  optimal  trajectories   nanating  from  a  single 
system  state  [B4]. 

The  solution  to  the  Riccati  equation  (8)  also  has  a  physical 
interpretation  which  will  be  useful  later. 

Lemma  2 . 2 

The  value  of  the  performance  index  for  the 
optimal  control  law  is 

J°(t    ,x    ,tj    =   |xTP(t    ,tjx  (10) 

oof  2o        ofo 

where  P  is  the  solution  to  the  matrix  Riccati 
equation  (8)  which  is  bounded  on  (t  , t  ) ,  i.e.^ 
there  exist:  no  conjugate  points. 

Proof 

t 


f 
J  =  2  /t  (xTQX  +  uTRu)dt 


2  Jt 
o 


Substitution  of  the  optimal  control  law, 


-1  T 
u  =  -R   G  Px, 


results  in 


J  =  |  J!   xT(Q  +  PGR  1GTP)x  dt.  (ID 


2 

Intersection  here  is  taken  so  as  to  exclude  the  case  of 

tangential  coincidence.   In  the  general  (non-linear)  case  extremals 
taken  sufficiently  close  (i.e.,  neighboring  trajectories)  must 
not  cross. 
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Define 

-1  T 
F,  =  F  -  GR   G  P 
k 

the  closed-loop  state  matrix  and  $,  (t  ,t)  as  its  corresponding  state 

k  o 

transition  matrix,  i.e. , 

$,  =  F,  §,   and   x(t)  =  $,  x  . 
k    k  k  k  o 

Equation  (11)  can  now  be  rewritten  as 

*f 

J  =  i  P   xT$J(Q  +  PGR~1GTP)§1  x  dt  .  (12) 

2  «t    o  k  k  o 

o 

Substitution  of  the  definition  of  F  into  the  Riccati  equation 

results  in 

-1  T 
P  =  -  F,  P  -  PF  -  PGR   G  P  -  Q 
k      k 


and  when  employed  in  (12)  to  replace  Q,  (12)  becomes 
J=zk       XX("PFk  "FkP  -^Vodt 


*f  --  *f 


J    =     -if.  XT§JW   X    dt     -    I     r  (XT$JP$,  X       +     XT§?W    X     )dt. 

2^t         okko  2^t  okko  okko 


Application  of  integration  by  parts  to  the  first  integral  above  results 
in  terms  which  cancel  the  second  integral,  leaving 

J  =  -  2  Xo5kPSkxo|t  =  I  XoP(*o'*f)Xo  "  I  xT(tf)P(tf,tf)x(tf). 


Since  P(tf,tf)  =  0,  this  is  the  desired  result. 

This  proof  is  considerably  different  from  the  standard  one 
(e.g. ,  [Al ,  p.  25])  in  that  it  does  not  require  the  use  of 


18 


Hamilton-Jacobi  theory.   This  lemma  not  only  provides  a  necessary  link 
in  the  solution  of  the  problem  considered  in  the  next  section  but 
provides  a  tie  with  dynamic  programming,  in  that  formulation  equation  (10) 
is  referred  to  as  an  "optimal  return  function"  [B5]. 

Sufficient  Conditions  -  Existence 

The  conjugate  point  condition  is  not  entirely  satisfactory  from 
the  standpoint  of  a  working  sufficiency  condition.   It  is  not  known  at 
the  onset  whether  or  not  the  optimization  problem  has  a  solution;  only 
after  the  complete  control  law  is  computed  is  one  assured  that  it  was 
not  all  for  naught.   What  is  required  then,  it  would  seem,  is  a  suffi- 
cient condition  somewhat  more  restrictive  than  the  conjugate  point 
condition  with  the  advantage  that  success  of  the  optimization  problem 
is  predetermined. 

Theorem  2. 3 

The  optimization  problem  of  minimizing  a  quadratic 
performance  index 


tf    T       T 
J  =  f    (x  Qx  +  u  Ru)dt,   R  =  R  positive  definite 


(13) 


subject  to  a  completely  controllable  linear  plant 

^   =  Fx  +  Gu  (14) 

dt 

always  has  the  global  minimum 

T  -1 

u(t)  =  -  K  (t)x(t)   where  K(t)  =  P(t,t  )GR 

and  ^£  _  -FTP  -PF  +  PGR_1GTP -Q,  P(t.,t.)=0  (15) 

dt  f   f 

for  all  finite  (t  ,t„)  if  Q  is  positive  semidefinite  and 

4.       ■  °  f 

symmetric. 
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Proof 

Suppose  that  Q  is  positive  semidef inite  and  the  Riccati  equa- 
tion (15)  diverges  at   t  =  t  <  t  .   A  solution  for  (15)  must  exist  in 
a  sufficiently  small  neighborhood  of  t  ,  further  P(t)  exists  for  all 
t  c(t+e,t)  ,  e>0;  hence,  by  Lemma  2.2 

J°(x(t1+e)  ,t  +e)  =  xT(tl+e)P(tl+e,tf)x(t1+e)  :>  0 ,        (16) 
for  e  <  t  -  t 

the  inequality  is  due  to  the  positive  definiteness  of  the  integrand 

of  J.   As  e  approaches  zero,  some  entry  of  P  becomes  unbounded.   It  can 

be  assumed  without  loss  of  generality  that  at  least  one  diagonal  element 

of  P  becomes  infinite;  otherwise,  some  2x2  principal  ninor  of  P  would 

be  negative,  contradicting  the  positive  semidef initeness  of  P  inferred 

in  (16).   Lete.  be  a  vector  which  is  all  zeros  except  for  the  ith  ele- 
i 

ment ,  a  one,  which  corresponds  to  a  diagonal  term  p..  of  P  which  becomes 
unbounded  as  €  approaches  zero;  then, 

J°(ei,t1+e)  =  Pii(t1+e,t  ). 

Since  J   is  the  optimal  performance  index,  a  performance  index 

resulting  from  any  arbitrary  control,  say  u  =  0,  must  be  greater. 

Then  if  t  +  e and e.   are  chosen  as  the  initial  time  and  state  of  the 

optimization  problem,  i.e. ,  t  =  t  +  e, 

o     1 

J°(e.  ,tl+e)=p.  .(tir€,tf)<  Jt   e^§T(T-to)Q5(T-to)e.  dT         (7) 


where   $(t-t  )  is  the  system  state  transition  matrix  and  $(t-t  )e. 
o  o   i 

is  the  resulting  free  trajectory.   Clearly,  the  integral  in  (17) 
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remains  bounded  over  the  finite  interval  (t  +e,t  ),  while  the  left- 
hand  side  becomes  unbounded  as  e  -*  0.   This  contradicts  the  original 
supposition  of  a  conjugate  point  at  t   and  proves  the  theorem. 

An  interesting  proof  of  this  theorem  using  Lyapunov's  second 
method  to  determine  the  stability  of  the  Riccati  equation  is  given  in 
[K2]. 

Throughout  the  literature  there  appear   sporadically  references 
to  theorems  such  as  2.3  as  necessary  and  sufficient  conditions  for  a 
solution  of  the  optimization  problem  to  exist  (e.g. ,  [D3,p.  557]); 
this  is  patently  false  and  numerous  counterexamples  exist  [G3]. 
No  necessary  and  sufficient  condition  for  the  non-existence  of  conjugate 
points  (other  than  integrating  the  Riccati  equation)  is  presently  known, 
although  many  researchers  are  actively  pursuing  these  conditions  and 
there  is  reason  to  believe  they  will  be  found  in  the  near  future  [B6]. 

2. 4  The  Optimal  Linear  Regulator  -  Infinite  Final  Time 

Existence  and  Stability 

The  Euler-Lagrange  equations  and  the  conjugate  point  condition 
apply  as  well  to  the  case  when  the  final  time  is  no  longer  finite. 
However,  as  before,  the  conjugate  point  condition  is  not  satisfactory 
as  an  existence  test;  the  following  theorem  extends  Theorem  2.3  to 
the  infinite  time  case. 
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Theorem  2. 4 

If  a  linear  plant  (14)  is  completely  controllable 
and  a  companion  quadratic  performance  index, 

j  =  f    (xTOx  +  uTRu)dt,   R  =  RT  positive  definite 

o 
has  Q  =  QT  positive  semidefinite  and  if  P(t,tf)  is  a 
solution  to  the  matrix  Riccati  equation  (15)  with 
P(tf,tf)  =  0,  then 

lira  P(t,tf)  -  P(t) 


(18) 


(19) 


exists  for  all  t  and  is  a  solution  of  (15). 

Proof  [K2] 

First,  it  must  be  shown  that  the  limit  (19)  exists  for  all  t. 
Since  the  plant  is  completely  controllable  for  every  XqJ  there  exists 
a  control  u'(t)  which  transfers  x  to  0  by  some  t  S  t^.   Set  u(t)  =  0 
for  t  >  tr   Then 


=  J(t  ,x  ,») 
o  o 


is  bounded  for  all  t>   tQ.   The  optimal  costs  are  also  non-decreasing 
as  t   -•.  Suppose  that  this  were  not  the  case,  that  is,  suppose 

^vvv^X-w-  for   t2>V 

where  u  (t)  is  the  optimal  control  corresponding  to  J°  and  u,,(t)  is 
the  optimal  control  resulting  in  J J.   Then  by  the  positive  definiteness 
of  the  integrand  of  (18),  use  of  control  u,,  will  result  in  a  lower  cost 
at  time  t   than  u  ,  thus  contradicting  the  optimality  of  l^  and 
demonstrating  the  assertion  that  the  optimal  costs  are  monotonic. 
The  limit  therefore  exists  for  all  t  by  the  well-known  result  that  all 
bounded  monotonic  sequences  possess  a  limit. 
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Now  it  is  necessary  to  demonstrate  that  (19)  is  a  solution  to 
the  Riccati  equation.   Define   P(t,t  ;A)  as  a  solution  of  (15)  with 
the  boundary  value  P(t   t  )  =  A.   Then,   using  the  continuity  of  solu- 
tions of  (15)  with  respect  to  its  boundary  values, 

P(t)    =      lim  P(t,t    ;0)    =      lim  P(t,t    •P(t1 ,t    -0)) 
t     -*c°  ^  t     ^co  112 

2  2 

=  P(t,t1;lim     P(t1,t2;0))    =   P(t, t    ;P(t   ) ) 
t2-co 

and  P(t)  is  a  solution  of  (15)  for  all  t. 

The  price  that  is  paid  for  guaranteeing  the  existence  of  a 
solution  over  all  time  (in  addition  to  requiring  Q  to  be  positive 
semidefinite)  is  the  complete  controllability  of  the  plant. 

Corollary  2. 1 

The  solution  P(t)  to  the  matrix  Riccati  equation  (15) 
for  the  optimal  linear  regulator  problem  of  Theorem  2.4  is 
constant  and  the  unique  positive  definite  solution  of 

PF  +  F  P  -  PGR~  GTP  +  Q  =  0,  (20) 

the  steady-state  (algebraic)  Riccati  equation. 

The  stationarity  of  P  follows  from  the  arbitrariness  of  the 

choice  of  the  initial  time  t  ,  since  all  initial  times  must  result  in 

o 

the  same  value  of  the  optimal  performance  index.   Equation  (20)  ensues 

when  the  effect  on  P  of  the  irrelevance  of  t   is  considered;  the  choice 

o 

of  the  positive  definite  solution  is  dictated  by  the  positive  definite- 
ness  of  (18)  and  uniqueness  is  guaranteed  by  the  conjugate  point  con- 
dition. 
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It  is  now  clear  that  the  optimal  linear  regulator  is  the  optimal 
control  theory  analog  of  the  classical  design  using  state  feedback. 

As  pointed  out  by  Kalman  [K2],  the  literature  contains  many 
references  where  the  stability  of  an  optimal  system  is  tacitly  assumed. 
For  example,  Letov  [Zl,  p.  378]  equates  (without  proof)  optimal  systems 
and  those  which  are  stable  in  the  sense  of  Lyapunov;  although  probably 
true  for  the  class  of  optimal  systems  which  are  also  stable,  it  is  not 
in  general.   This  can  be  easily  demonstrated  with  a  simple  example. 
Let  a  scalar  input  plant, 

x  =  Fx  +  gu, 
have  one  or  more  eigenvalues  with  positive  real  parts  and  let  the  per- 
formance index  to  be  minimized  be, 


J  =  /t   u  dt, 


that  is,  Q  =  0  and  R  =  1.   The  problem  so  defined  clearly  has  a  solu- 
tion (Theorem  2.3),  which  is 

u(t)  =  0. 
Stability  was  not  a  result  of  optimization  in  this  case  primarily 
because  the  states  are  not  reflected  in  the  cost  (performance  index) ; 
had  the  plant  been  stable,  however,  the  optimal  system  would  have  been 
also.   The  final  theorem  in  this  section  formalizes  this  observation 
into  a  general  result. 
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Theorem  2.5  [Kl] 

An  optimal  linear  regulator  problem  satisfying  the 
conditions  of  Theorem  2.3  results  in  an  asymptotically 
stable  control  law  if 

T 
i)  the  pair  [F,H],  where  Q  =  HH   is  completely 

observable,  i.e.,  [H,FTH, . . . (FT)n-1H]  has 
full  rank; 

and  only  if 

ii)  the  linear  subspace 

X  =  {x*  0|  ||HeFtx||=  0} 

Ft 
of  the  state  space  is  null  (i)  or  e  x,  x  e  X 

is  an  asymptotically  stable  response  in  the 
sense  of  Lyapunov. 


Proof 
i)   By  the  assumption  of  observability  and  the  positive  semi- 
definiteness  of  Q,  the  integrand  of  the  performance  index  must  be  posi- 
tive along  any  non-zero  trajectory  of  the  plant.   Then 

J°(t,x(t))  =  xT(t)Px(t)  >  0,   x(t)  ^   0  (21) 

along  optimal  trajectories.   Differentiation  of  the  performance  index, 

00 

J(t,x(t))  =  f  (xTQx  +  uTRu)dT, 

T       T 
results  in     J(t,x(t))  =  -  (x  Qx  +  u  Ru)<  0,  x(t)  ^  0,  which  is 

negative  along  all  non-trivial  trajectories.   Then  (21)  is  a  Lyapunov 

function  [L2]  and  the  system  is  asymptotically  stable  by  Lyapunov' s 

Second  Stability  Theorem  [L3,  p.  37], 

ii)   When  certain  of  the  states  are  not  observable  in  the  cost, 

Lyapunov1 s  method  is  effectively  applied  only  to  a  subsystem  (i.e. ,  the 
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states  observed  in  the  cost)  and  the  more  general  result  follows  when 
the  stability  of  the  remainder  of  the  system  is  considered. 

The  requirements  for  stability  placed  on  Q  by  Theorem  2.5  are 
not  unduly  restrictive  in  that  they  reflect  good  engineering  judgment; 
that  is,  if  a  state  has  significant  effect  on  system  performance  it 
should  influence  the  design  process.   Some  authors  choose  to  include 
stability  tacitly  in  the  definition  of  the  performance  index  by  con- 
sidering only  the  case  where  Q  is  positive  definite  (e.g. ,  [SI]). 
Stability  is  assured  since  a  positive  definite  Q  must  have  a  non- 
singular  factor  and  condition  (i)  of  Theorem  2.5  is  clearly  satisfied. 

2 .  5   Summary 

Over  the  past  ten  years  the  optimal  linear  regulator  problem 
has  become  one  of  the  most  widely  studied  optimal  control  problems. 
Perhaps  the  best  way  to  summarize  this  unique  problem  is  to  review  the 
qualities  which  set  it  apart. 

1.  The  solution  to  the  problem  is  always  approached  in  the 
same  fashion  regardless  of  the  specific  system  or  weight- 
ing matrices  employed  and  in  contrast  to  the  majority  of 
optimal  control  problems,  the  solution  is  explicit. 

2.  In  a  restricted,  but  very  large,  set  of  weighting  matrices 
(i.e. ,  Q  positive  semidef inite)  the  solution  is  guaranteed 
to  exist. 

3.  The  solution  is  always  in  the  form  of  a  linear,  constant 
feedback  control  law,  as  opposed  to  the  general  optimal 
control  problem  where  a  feedback  formulation  is  not  directly 


26 


obtainable.   With  another  minor  concession  to  generality 
(Theorem  2.5)  asymptotic  stability  will  be  assured  from 
the  onset. 
4.   Numerical  computation  of  the  control  laws  is  both  straight- 
forward and  relatively  inexpensive  (in  terms  of  computing 
time)[B7];  in  addition,  a  multitude  of  numerical  schemes 
are  available  [Bl]. 
Under  suitable  assumptions  of  continuity  and  corresponding 
definitions  of  controllability  and  observability  [D2]  most  of  these 
remarks  are  also  applicable  to  linear  time-varying  systems  and  non- 
constant  weighting  matrices  [K2]. 


CHAPTER  III 
THE  INVERSE  OPTIMAL  LINEAR  REGULATOR  PROBLEM 

3. 1   Introduction 

Optimality  in  itself  is  not  necessarily  a  desirable  character- 
istic for  a  system  to  possess.   For  instance,  suppose  it  is  desired 
to  minimize  the  sensitivity  of  an  existing  system's  response  to  com- 
ponent variations  by  use  of  an  appropriate  controller.   A  system 
design  which  obviously  has  a  minimal  sensitivity  is  one  that  does 
nothing  whatever;  then  one  candidate  for  an  optimal  controller  would 
be  one  that  turns  the  system  off.   In  this  case,  the  optimal  design  would 
probably  not  be  satisfactory  because  it  was  not  implicit  in  the  procedure  that 
the  resulting  system  should  perform  in  some  acceptable  manner  in  addi- 
tion to  minimizing  the  performance  index. 

The  lesson  is  clear,  optimality  may  be  a  frivolous  notion 
unless  its  ramifications  are  thoroughly  understood  in  the  context  of 
total  system  behavior.   In  the  present  chapter  the  implications  of 
optimality  are  considered  for  a  single  input  linear  plant  subject  to 
a  quadratic  performance  index  to  set  the  stage  for  its  use  as  a  viable 
design  tool. 

The  system  to  be  considered  in  this  chapter  can  be  described 
by  the  state  equation, 

fg=Fx+gu,  (1) 

dt 
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and  occasionally  it  will  require  an  associated  feedback  control  law, 

u  =  -  kTx,  (2) 

where  F  is  a  constant  n  x  n  matrix  and  g,k  are  constant  n-vectors. 
In  conjunction  with  this  plant,  a  quadratic  performance  index, 

CD 

J*   T        2 
(x  Qx  +  ru  )dt,  (3) 

o 

will  be  studied.   The  control  weighting  factor,  r,  will  be  taken  with- 
out loss  of  generality  as  unity  throughout  the  remainder  of  this  work. 

3.2   When  Is  a  Linear  Control  System  Optimal? 

The  original  impetus  to  investigate  the  "inverse  optimal  control 
problem"  is  generally  considered  to  be  a  1964  paper  by  R,  E.  Kalman, 
bearing  the  same  title  as  this  section  [Kl],  although  it  has  a  much 
older  history  in  the  Calculus  of  Variations.   In  that  paper,  Kalman 
considers  the  theoretical  criterion  for  a  linear  system  to  be  optimal 
with  respect  to  a  restricted  class  of  quadratic  performance  indices. 
His  principal  result  is  reviewed  next 

Theorem  3. 1  [Kl] 

Consider  a  completely  controllable  linear  plant 
(1)  with  a  stable  completely  observable  control  law  k 
(2).   Then  k  is  an  optimal  control  law  with  respect 
to  a  quadratic  performance  index  (3) ,  with  Q  positive 
semidef inite ,  if  and  only   if, 


|l   +    kT§(juo)g|2    S   1  (4) 


for    all    real   tu ,    where      $(s)    -    (si    -   F) 
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Proof 
i)   Necessity.     If  k  is  an  optimal  control  law,  then  by 
Theorems  2.3  and  2.4 

k  =  Pg  (5) 

where  P  is  a  solution  to  the  steady-state  Riccati  equation, 

PF  +  FTP  -  PggTP  +  Q  =  0.  (6) 

Substituting  (5)  into  (6)  permits  it  to  be  rewritten  as 

PF  +  FTP  -  kkT  +  Q  =  0.  (7) 

Adding  and  subtracting  sP  to  (7)  results  in 

-P(sl  -  F)  -  (-si  -  FT)P  -  kk  +  Q  =  0. 

T  T 
Premultiplying  and  postmultiplying  the  above  by  g  §  (-s)  and  $(s)g, 

respectively  (where  $(s)  is  as  defined  in  equation  (4))  and  substitut- 
ing (5)  where  appropriate,  leads  to 

g  $  (-s)k  +  k  §(s)g  +  g  $  (-s)kk  $(s)g  =  g  5  (-s)Q$(s)g 
which  may  be  factored  as 

[1  +  gT$T(-s)k][l  +  kT$(s)g]  =  1  +  gT$T(-s)Q$(s)g.  (8) 

Noting  the  semidef initeness  of  Q  proves  the  necessary  condition.  . 

ii)  Suf f iciency  Proof  of  the  sufficient  part  of  the  theorem 
is  not  particularly  enlightening  and  requires  concepts  that  are  yet 
to  appear,  hence  it  will  be  deferred  to  the  appendix. 

3.3   Implications  of  Optimality 

It  is  somewhat  surprising  that  the  condition  for  optimality  (4) 
should  be  a  frequency  domain  criteria  in  that  all  of  the  previous 
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analysis  of  the  problem  (Chapter  II)  was  carried  out  in  the  time  domain. 
The  study  of  optimal  systems  in  the  frequency  domain  may  be  responsible 
for  the  current  interest  shown  by  Kalman  I  KG  1  and  others  in  realization 

e  >ry  from  the  point  of  view  of  invariants  and  thus  promises  to  have 
far  greater  i  nfluence  in  years  to  come. 

Examination  of  the  optimality  criterion  (4)  reveals  that  the 
quantity  whose  modulus  is  taken  on  the  left-hand  side  of  the  equation, 
i.e., 

fk(ju>)  =  1  +  k  §(ju))g 

is  the  "return  difference"  of  classical  feedback  theory  [B8].   It  can 
be  interpreted  as  the  difference  between  a  unity  input  (in  the  frequency 
domain)  to  the  system  and  what  would  consequently  be  returned  as  feed- 
back.  Graphically,  this  is  illustrated  in  Figure  3.1  as  the  difference 
between  the  input  at  node  a  and  the  feedback  at  node  b. 


Plant 


a 
0— 

1- 


$(jw)g 


Feedback 
Gains 


f,  =  a  -  b 
k 


Figure  3. 1   Return  Difference 


The  requirement  of  (4)  that  the  modulus  of  the  return  differ- 
ence be  greater  than  unity  is  the  celebrated  result  of  classical 
feedback  theory  that  the  sensitivity  of  system  response  to  variations 


in  plant  parameters  is  reduced  by  the  addition  of  feedback  [B8] . 
Further  analysis  of  inequality  (4)  reveals  that  it  requires  that  the 
Nyquist  locus  avoid  a  circle  of  unit  radius  about  the  point  (-1,0) . 
Figure  3.2  pictures  a  Nyquist  plot  of  a  hypothetical  optimal  system 
constructed  from  a  plant  with  two  eigenvalues  with  positive  real  parts. 


Figure  3.2    Nyquist  Plot  of  Optimal  System 


Two  additional  observations  may  be  immediately  made  from 
Figure  3.2.   First,  the  phase  margin  of  an  optimal  system  must  be  at 
least  60°,  since  the  closest  points  (in  phase)  to  the  negative  real 
axis  that  the  locus  can  cross  the  unity  magnitude  contour  (points  A 
and  A'  on  Figure  3.2)  are  displaced  60   from  the  axis.   Second,  the 
system  will  clearly  remain  stable  regardless  of  how  the  feedback  gain 
is  increased.   Then  if  gain  margin  is  defined  as  that  factor  by  which 
feedback  gain  may  be  increased  before  system  instability  occurs  [E2] , 
the  gain  margin  of  an  optimal  system  may  be  said  to  be  infinite. 
However,  if  gain  margin  is  defined  as  the  reciprocal  of  the  magnitude 
of  the  loop  gain  at  phase  crossover,  i.e., 
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r„   =   i where      Arg    rkT$(jvu    )g]    =    180°, 

i  T  I  C 

|k  *(ju>c)g| 

then  an  optimal  system  has  a  gain  margin  of  at  least  -  or  approximately 

-6  db.   Kuo  indicates  that  the  two  definitions  given  are  equivalent 

[K7,  p.  398],  but  consideration  of  the  example  of  Figure  3.2  reveals 

that  this  is  not  the  case  when  the  plant  is  non-minimum  phase. 

A  third  definition,  and  the  one  that  should  be  understood  in 

any  references  to  gain  margin  in  the  sequel,  is  the  following  fT2]  : 

Let  g.  be  the  factor  by  which  the  feedback  gain  may  be 
increased  until  instability  occurs  and  g   be  the  factor 
by  which  the  gain  may  be  decreased  and  stability  main- 
tained, then 


GM 


smaller  of  fg.,l/g,] 


This  definition  is  more  reasonable  in  that  it  reflects  the  actual 
gain  disturbance  required  for  instability  and  is  easily  determined 
from  the  Nyquist  diagram.   The  gain  margin  for  an  optimal  regulator, 
in  the  sense  of  Theorem  3.1,  then  is  at  least  6  db. 

It  is  instructive  at  this  point  to  compare  these  figures  of 
gain  and  phase  margin  with  those  required  in  practice.   Jones,  Moore 
and  Tecosky  in  Truxal's  classic  handbook  [T2.PP-  19-14]  recommend 
phase  margins  of  40°  to  60°  and  gain  margins  of  at  least  10  db   for 
applications  typical  of  chemical  process  control.   Generalities  of 
this  sort  are  of  course  subject  to  severe  criticism  but  nonetheless 
it  can  be  concluded  that  minimum  optimal  regulator  specifications 
compare  favorably  with  those  desired  of  classical  designs. 
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3. 4   Companion  Matrix  Canonical  Form 

It  is  well  known  that  for  any  non-de.~ogacory  matrix,  F,  there  exists 
a  similarity  transformation,  T,  such  that 

-1 

F  =  TFT 

is  the  companion  matrix  [Ml]  of  F.   A  companion  matrix  being  a  matrix 

of  the  form: 

0     1    0  ...  0 
0     0     1  ...  0 


F  = 


(9) 


where 


Lai  "a2  -V-'-Vj 


n  k  1 

V(s)  =  s  +   £  as    =  det  (si  -  F>. 

k=l  k 


That  is,  the  last  row  of  the  companion  matrix  of  F  contains  the  nega- 
tive of  the  normalized  coefficients  of  the  characteristic  polynomial 
of  F  (and  F) ,  with  the  remainder  of  the  matrix  null  save  'a  superdiag- 
onal  of  ones.  Some  authors  choose  to  refer  to  the  transpose  of  the 
matrix  defined  above  as  the  companion  matrix  [G4]  jn  the  following 
either  definition -will  suffice;  however,  (9)  will  be  assumed  for  con- 
sistency. 

When  F  assumes  the  role  of  the  state  matrix  in  a  completely 
controllable  linear  system  (1)  the  similarity  transformation,  T,  which 
transforms  F  to  its  companion  matrix,  places  the  system  as  a  whole  in 
a  canonical  form.   Specifically,  if  z  =  Tx  represents  the  change  of 
basis  described  above  and  the  linear  completely  controllable  system  is 


:;■! 


x  -  Fx  +  gu 
T 

y  =  II  x, 


then  |K1] 


z  =  Fz  +  gu 

y  =  H  z 

F  =  TFT   ,  the  companion  matrix  of  F 

0 
0 

g  =  Tg=   :   ,  and  H  =  (T-1)TH. 

0 
1 

Complete  controllability  allows  the  specified  form  of  g  to  be  chosen. 
Since  there  is  no  loss  of  generality,  controllable  systems  will  be 
assumed  in  this  coordinate  system  and  the  "hat"  notation  will  be  omitted. 

The  problem  of  actually  computing  the  required  transformation 
has  lately  occupied  a  good  deal  of  the  literature,  e.g.  ,  [~J1,R2,R3]. 
Computation  of  the  transformation  is  conceptually  not  very  difficult 
and  some  straightforward  numerical  solutions  have  been  formulated 
Tesp.  R3] . 

The  companion  matrix  canonical  form  is  often  referred  to  as  the 
phase-variable  canonical  form  fS2];  this  cognomen  alludes  to  the  prop- 
erty that  each  state  is  the  derivative  of  the  preceding  state. 
A  second  property  that  will  find  application  is  the  fact  that  g  lias 
only  one  non-zero  element,  hence,  the  control  directly  affects  only 
the  last  state  (highest  derivative). 
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3. 5   Characterization  of  the  Equivalence 
CI  ass  of  Q' s 

Employment  of  the  coordinate  system  reviewed  in  the  last  section 

permits  additional  insight  to  be  gleaned  from  Theorem  3.1.   First  it 

will  be  necessary  to  exploit  a  characteristic  of  the  companion  matrix 

canonical  form. 

Lemma  3. 1 

For  a  completely  controllable  linear  system  in  the 
companion  matrix  canonical  form  (1) 


S(s)  =  Cp(s)8(s)g  = 


.n-1 


where   $(s)  =  (si  -  F) 

and    cp(s)  =  det  (si  -  F)  . 

Proof 

To  prove  the  theorem  it  is  sufficient  to  show  that 

(sI-F)cp(s)$(s)g  =  (sI-F)S(s)  =  Cp(s)g. 
With  F  in  companion  matrix  form 


(sI-F)S(s) 


-1    0 
s   -1 


,n-l 


su+.E,a.s 
i=l  l 


i-1 


which  is  clearly  cp(s)g  and  the  theorem  is  proven. 

This  lemma  allows  the  simplification  of  equation  (8)  used 


the  proof  of  Theorem  3. 1  to 


1  + 


kTS(juj) 
9(juj) 


1  + 


S  (-ju-)QS(jtu) 
|<P(Jw)|2 


(10) 
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Multiplying  by  |cp(j(ju)  |   results  in 

|cp(juu)  +  kTS(joj)|   =  |cp(ju))|   +  ST(-jtu)QS(juj). 

With  the  observation  that  the  expression  inside  the  absolute  value  marks 
is  the  closed-loop  characteristic  polynomial  (which  will  be  denoted 
cp  (s))  permits  further  simplification  to 

|cpkOuO|   -  |cp(jcu)|   =  S  (-jio)QS(ju)).  (11) 

This  relation  (11)  defines  a  polynomial  by  which  the  open-  and 
closed-loop  characteristic  equations  of  an  optimal  system  must  be 
related  to  the  state  weighting  matrix  of  a  corresponding  performance 
index.   This  polynomial  will  appear  frequently;  consequently,  it  will 
be  convenient  to  refer  to  it  with  the  following  notation: 

Y(u>)  =  |9k(jw)|2  -  fcp  Cjou)  | 2  =  ST(-juj)QS(juu),  (12) 

or  occasionally  as 

Y(Q;u>)  =  ST(-j(ju)QS(juu) 

to  emphasize  the  functional  relationship  of  Q. 

Recalling  that  characteristic  polynomials  are  invariant  under 
similarity  transformation  and  the  assumption  of  positive  semidef initeness 
for  Q  leads  to  a  useful  corollary  of  Theorem  3. 1. 

Corollary  3.  1 

A  completely  controllable  scalar  linear  system 

(1)  with  a  completely  observable  stable  control  law  k 

is  optimal  with  respect  to  a  quadratic  performance  index 

(2)  with  Q  positive  semidefinite  and  completely  observ- 
able if  and  only  if 

|cp  (juj)  |   -  |tp(j'Jj)|   ^  0   for  all  real  u), 
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where  cp(s)  =  det  (sI-F)  and  cp  (s)  =  det  (si  -  F  +  gk  ), 
the  open-  and  closed-loop  characteristic  polynomials, 
respectively. 

The  requirement  that  the  feedback  control  be  completely  observ- 
able again  appears  without  apparent  justification.   Suppose  that  k  were 
not  completely  observable,  then  the  pair  of  polynomials  which  comprise 
the  expression 


.T,,.  ,     K(ju)) 

k  $(ju))g  = 


will  have  a  common  factor  IK5].   This  can  be  thought  of,  in  the  clas- 
sical control  sense,  as  having   a  zero  on  top  of  a  pole  which  prevents 
the  response  of  the  pole  from  being  observed  at  the  output.   Let  the 
common  factor  be  yU^)  ancl  denote  the  polynomials  with  y   factored  with 

a  prime,  then  equation  (10)  becomes 

2 


K  (jcu)y(juj) 


cp  (juj)y(Ju)) 


=  1  + 


S  (-jiu)QS(jcu) 


Jcp'CjU))  +K/(ju))|  |y(jiu)|   =  |cp'(Jw)|  |yUu>)|  +  ST(-ja))QS(juj). 

T 
The  implication  is  that  S  (-jw)  QS(joj)  must  contain  the  common  factor 

Iy(Jw)|   as  well.   Since  Q  is  positive  semidef inite,  it  can  be  factored 

T 
as  Q  =  HH  ;  then  by  the  preceding  statement  the  vector  of  polynomials 

formed  from 


Through  a  minor  abuse  of  the  language,  the  terms  "completely 
observable  Q"  and  "observability  from  the  cost"  should  be  taken  to 
mean  that  the  pair  [F  H]  are  completely  observable  when  H  is  any  factor 
of  Q  such  that  Q  =  HH.   Similarly,  "completely  observable  k"  means 
that  the  pair  [F,k]  are  completely  observable. 
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H  §(ja>)g 


(13) 


may  have  the  common  factor  y(Ja))  and  the  Q  may  not  be  completely  observ- 
able.  The  reason  for  the  uncertainty  in  the  observability  of  Q  is  the 
ambiguity  in  the  location  of  the  zeros  of  the  factor  of  (13)  correspond- 
ing to  |Y(jtu)  I  .   Depending  on  the  choice  of  H,  (13)  can  have  zeros  in 
the  right  or  left  half-plane  which  may  or  may  not  eclipse  a  pole  of  the 
system.   If  k  is  completely  observable,  however,  there  will  be  no  common 
factor  regardless  of  how  H  is  chosen  and  any  Q  which  satisfies  (13)  will 
be  completely  observable.   This  is  best  illustrated  with  an  example. 

Example  3. 1 

Consider  the  second  order  plant, 


F  = 


0    1' 


-1 


with  the  feedback  law, 


The  feedback  is  clearly  not  observable  and  in  fact  cancels  one  of  the 
plant  poles  at  -1  as  shown  in  Figure  3.3. 


-2 


Iiu 


x   -   plant   poles:  -1,    -1 

D  -   closed-loop   poles:    -1,    -2 


Re  - 


Figure  3.3   System  Pole  Locations 
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The  control  law  is  nonetheless  optimal  and  three  performance  indices 

.   2 
which  are  minimized  by  it  have  Q  s, 


and 


Ql  = 

"3      3^ 
_  3      3-j 

hx   =  J5 

1 
-  1- 

Q2   = 

"~  3    -3" 
_-3      3_ 

s-jQ_ 

-3   <r 

Q3   = 

^-0      3_ 

The  first  is  not  observable;  (13)  has  a  zero  at  -1  which 
coincides  with  a  plant  pole.   The  second  is  identical  to  Q  except 
that  the  zero  of  (13)  is  reflected  about  the  axis  (at  +  1)  and  is 
consequently  observable.   The  last  is  non-singular  and  avoids  the 
difficulty  entirely.   The  important  point  to  note  is  that  the  control 
is  optimal  for  all  of  these  Q's.   Complete  observability  from  the  cost 
is  not  required  for  optimality;  it  is  only  necessary  that  any  plant 
poles  which  are  unobservable  appear  in  the  closed-loop  system  as  well. 
That  is,  optimal  feedback  cannot  move  system  poles  which  are  unobserv- 
able in  the  cost. 

The  requirement  that  k  be  completely  observable  is,  in  a  sense, 
an  "inverse"  sufficient  condition  to  the  observability  of  Q;  if  the 
condition  is  met,  then  any  Q  for  which  the  system  is  optimal  must  be 
completely  observable.   If,  however,  k  is  not  completely  observable  but 
Y  is  non-negative,  there  may  still  exist  one  or  more  positive  semidef i- 
nite  Q's  which  are  completely  observable  for  which  the  system  is  optimal, 


These  statements  will  be  justified  later. 
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Another  corollary  to  Theorem  3. 1  results  from  consideration  of  the 
preceding  remarks. 

Corollary  3.2 

A  completely  controllable  scalar  linear  system  with 
a  stable  feedback  control  law  is  optimal  with  respect  to 
a  quadratic  performance  index  with  Q  positive  semidefinite  if 

|9  (juu)  f   -  |cp(jcju)|   £  0     for  all  real  u), 


and 


a)  cp(s)  and  Cp  (s)  are  relatively  prime  which 
further  insures  that  any  Q  for  which  the 


control  law  is  optimal  is  completely  observable, 

or 

b)  the  plant  is  stable. 

Physically,  the  (a)  condition  can  be  interpreted  as  meaning  that 
the  aforementioned  observability  difficulties  do  not  arise  if  the  feed- 
back moves  all  the  plant  poles.   This  is  the  case  because  the  common 
factor  in  the  discussion  preceding  the  corollary  does  not  exist.   The 
second  (b)  condition  removes  the  ambiguity  by  restricting  the  plant 
poles  to  the  left-half  complex  plane  (including  the  imaginary  axis) 
and  no  difficulties  of  the  sort  discussed  will  occur.   Some  necessary 
and  sufficient  conditions  similar  to  those  of  Corollary  3.2  will  be 
considered  coincidentally  with  later  results. 

In  order  to  provide  greater  insight  into  the  optimal ity  condi- 
tion in  a  strictly  mathematical  sense,  it  will  be  convenient  to  consider 
the  requirements  for  Q  implied  by  equation  (11).   In  the  same  fashion 
that  the  two  preceding  corollaries  dealt  with  the  left-hand  side  of  the 
equation  and  its  relation  to  the  system,  it  is  desirable  to  consider 

what  the  right-hand  side,  i.e. , 

Y(Q;uj)  =  ST(-ju))QS(juj)  , 
portends  for  the  corresponding  performance  index. 
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Before  proceeding,  two  well-known  lemmas  on  factorization  will 
be  stated;  the  first  concerns  a  factorization  of  real,  even  polynomials 
and  the  second  a  factorization  of  positive  semidefinite  matrices. 

Lemma  3.2 

If  and  only  if  T(oj)    is  a  real,  even  polynomial  and 

r(uu)  ^  0    for  all  real  U), 

then  there  exists  a  unique  "spectral"  factor  y(s)  such 
that 

and  y(s)  has  only  zeros  with  non-positive  real  parts. 

Many  proofs  of  this  lemma  are  recorded  in  the  literature;  for 
example,  see  Brockett  [B2,  p.  173  f f ] . 

Lemma  3.3 

A  real,  symmetric  matrix  Q  may  be  factored  as 

T 
Q  =  HH  , 

where  H  is  a  real  matrix  of  rank  (H)  =  rank  (Q)  if 

and  only  if  Q  is  positive  semidefinite. 

This  is  a  fundamental  theorem  which  is  encountered  in  the  study 
of  quadratic  forms;  see  [A2 ,  p.  139]. 

The  groundwork  has  now  been  laid  to  present  a  theorem  which 
concisely  places  the  system-theoretic  results  of  Theorem  3. 1  into  a 
mathematical  framework. 

Theorem  3. 2 

For  every  real,  even  polynomial  T(u))  of  order 
2(n-l),  there  exists  a  real,  symmetric,  positive  semi- 
definite,  nth  order  matrix  Q  such  that 

Y(Q;uj)  =  r(cu)  , 

(Y  as  defined  in  (12))  if  and  only  if, 

F(uj)  ^  0     for  all  real  uu. 
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Proof 
i)   Sufficiency:   If  r(uu)  is  real,  even  and  non-negative, 
Lemma  3.2  assures  that  a  real  factor  y(s)  exists  such  that 

Y(Jw)y(-Jw)~=  r(cu).  (14) 

Let  h  be  a  real  vector  composed  of  the  coefficients  of  y(s)  ordered 
with  the  constant  term  first  and  the  coefficient  of  the  (n-l)st  term 

last.   Then,  clearly, 

T  T 

Y(Ju>)  =  h  S(jtu)   and  y(-j(Xi)    =    S  (-juu)h, 

where  S  is  as  defined  in  equation  (12)  of  Lemma  3.2.   The  product  (14) 

can  be  identified  as 

Y(hhT;ai)  =  ST(-juu)hh  S(jou)  =  Y<-Ja)>Y(Jcu)  =  r(tu) 

T 
and  hh   provides  a  positive  semidefinite  Q  which  satisfies  the  suffi- 
ciency part  of  the  theorem. 

ii)   Necessity:   It  must  be  shown  that  any  positive  semidefi- 
nite Q  results  in  a  Y(Q;uo)  that,  is  a  real,  even,  non-negative  poly- 
nomial.  Since  Q  is  positive  semidefinite,  the  Hermitian  form, 

T(u))  =  Y(Q;u))  =  ST(-joj)QS(jco), 
is  non-negative  and  is  clearly  real  and  even.   The  theorem  is  proved. 

Theorem  3.2  relates  the  optimality  condition  of  Corollary  3.1  to 
a  realizability  condition  for  positive  semidefinite  Q' s.   Corollary  3.1 
states  that  if  the  polynomial 

Y(uj)  =  |cpk(jw)|   -  jcpCJcu)  | 
is  non-negative  and  the  feedback  control  law  is  stable,  the  closed-loop 
system  is  optimal.   Theorem  3.2  states  that  if  an  arbitrary  real  even 
polynomial  is  non-negative,  then  there  exists  a  positive  semidefinite 
matrix  Q  that  generates  the  polynomial  by  equation  (12).   If  the  system 
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is  in  companion  matrix  canonical  form,  the  polynomial  studied  for 
optimal ity  and  the  polynomial  of  Theorem  3.2  are  one  and  the  same  as 
revealed  by  equation  (11).   Together  they  constitute  all  that  is 
required  to  determine  the  optimality  of  a  given  closed-loop  system 
configuration  with  respect  to  a  positive  semidefinite  Q  and  to 
construct  a  corresponding  performance  index. 

The  proof  of  Theorem  3.2  provides  a  "recipe"  for  computing 
at  least  one  performance  index  which  is  minimized  by  a  given  canonical 
optimal  system;  that  is,  by  spectral  factorization.   This  realization 
also  prevents  the  performance  index  from  obscuring  unstable  system 
poles  and  thus  avoids  the  observability  problem.   At  this  point,  it 
would  be  instructive  to  return  to  Example  3. 1  and  verify  system 
optimality  and  the  choices  of  Q. 

Example  3. 2 

Consider  the  second  order  plant  and  feedback  law  of  Example  3.1 


F  = 


0    1 
.-1   -2 


k  = 


The  open-  and  closed-loop  characteristic  polynomials  are 

2 

tp(s)  =  s  +  2s  +  1, 

2 
cp,  (s)  =  s   +  3s  +  2, 
k 


and 


Y(u>)  =  |cp  (ju))|   -  |cp(ju))|   =  3ju2  +  3, 


is   clearly  non-negative.   Then  since  the  plant  is  stable,  the  opti- 
mality of  the  feedback  law  follows  directly  from  Corollary  3.2b. 
Theorem  3.2  insures  that  a  positive  semidefinite  Q  can  be  found  which 
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will  form  a  performance  index  minimized  by  the  control  law, and  the 
proof  hints  at  how  one  such  Q  can  be  constructed.   Following  the  proof 
the  spectral  factor  of  Y  is  computed,  i.e.  , 

i|f(s)  =  ^3(s+l)   and 

Y  =  <r(-ju))i|f(jaj). 
Now  a  vector  h  is  constructed  composed  of  the  coefficients  of  l|r(s) 
ordered  with  the  constant  term  first, 


h  =  v^ 


and   Q  =  hh 


"3   3-| 
_3   3 


which  is  the  choice  for  Q   in  Example  3. 1. 

In  order  to  avoid  an  increasing  awkwardness  in  notation,  the 
following  definition  is  required. 

Definition  3. 1 

A  state  weighting  matrix  Q„  is  said  to  be  equivalent 
to  another  weighting  matrix  Q^  if  for  a  given  linear  system  (1) 
the  quadratic  performance  indices  (3)  formed  from  Qi  and  Qo 
are  minimized  by  the  same  control  law.   This  equivalence 
will  be  denoted  by  a  tilde,  i.e. , 

Ql~  %• 
and  the  set  of  all  matrices  which  are  equivalent  for  a 
given  system  and  optimal  control  law  will  be  referred  to  as 
the  "equivalence  class  of  Q' s"  for  that  optimal  system. 


Note  that  this  equivalence  relation  is  correctly  defined  in  the 
strictest  sense  [A3],  that  is,  it  is: 
i)   reflective   -  Q  ~  Q 
ii)   symmetric    -  if  Q  ~  Q.,  then   Q  ~  Q 
iii)   transitive   -  if  Q  ~  Q.  and  Q  ~  Q  ,  then  Q  ~  Q 
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The  definition  fails  to  draw  a  distinction  between  symmetric  and 
non-symmetric  matrices.   Consistent  with  the  remainder  of  this  work, 
symmetric  weighting  matrices  will  be  tacitly  assumed.   There  is  no  loss 
of  generality  in  the  assumption  of  Q  symmetric;  if  A  is  a  non-symmetric 


matrix,  then  it  is  easy  to  see  that  it  can  be  replaced  in  a  quadratic 

1 


1       T 
form  by  the  symmetric  matrix  —(A  +  A  )  without  altering  the  value  of 


the  form  I"H1] . 

The  next  theorem,  which  may  be  considered  a  central  result  of 
this  section,  prescribes  how  all  quadratic  performance  indices  minimized 
by  a  given  system  are  related. 

Theorem  3. 3 

For  a  scalar  system,  x  =  Fx  +  gu,  in  companion  matrix 
canonical  form  with  a  stable  optimal  control  law  k,  a  state 
weighting  matrix  Q2  is  equivalent  to  a  weighting  matrix  Q^ 

which  forms  a  quadratic  performance  index  minimized  by  k 
if  and  only  if 

a)  YCQgju))  =  Y(Ql;u;), 

b)  the  linear  subspace  of  the  state  space, 

T 
T  F  t   Ft 


X  =  fx  4   0  |x  e    Q2  x  =  0} 


Ft 
is  null  or  e   x  for  x  e  X   is  an  asymptotically  stable 

response  in  the  sense  of  Lyapunov,  and 

c)   the  optimization  problem  for  ^   and  the  given  system 
has  no  conjugate  points  on  t  e  [0,<=°). 

Proof 
i)   Sufficiency:   By  Theorem  2.1  and  (c)  above,  the  minimiza- 
tion of  the  performance  index  with  Q  results  in  a  unique  optimal  control 
law  which  by  Theorem  2.5  and  (b)  is  also  stable;  denote  this  control  law 
as  k2.   Suppose  that  k2  /  k   where  k   is  the  optimal  control  law 
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corresponding  to  Q  .   Then,  since  the  coefficients  of  the  closed-loop 

characteristic  polynomials  are  the  sum  of  the  plant  characteristic 

polynomial  coefficients  and  the  entries  of  k   (or  k  ) , 

cp   (s)  4   cp   (s)  , 
k2       kl 

and  by    the   stability   of   k      and  k 
'•1.2 

|cp      (ju))|2  *    |cp      (joj)|2 
K2  kl 

(i.e.,  both  characteristic  equations  have  zeros  in  the  left  half-plane). 

This  last  inequality  with  the  definition  of  Y  (12)  contradicts  part  (a) 

of  the  hypothesis  and  sufficiency  is  demonstrated. 

ii)   Necessity:   It  must  now  be  shown  that  if  Q  ~  Q  ,  then 
(a) ,  (b) ,  and  (c)  follow.   Parts  (b)  and  (c)  are  implicit  in  that  they 
are  necessary  and  sufficient  conditions  for  stability  and  existence, 
respectively,  of  the  optimization  problem  with  Q  .   The  first  part 
results  from  the  equality  of  k   and  k   and  the  definition  of  Y  (12). 

This  theorem  could  easily  have  been  rewritten  to  include  the 
case  where  the  resulting  control  law  was  not  stable;  however,  this 
complication  would  have  no  usefulness  and  would  obscure  insight 
into  the  mechanism  of  equivalent  Q' s.   Many  of  the  results  to  follow 
can  be  extended  to  include  unstable  optimal  control  laws  but  any 
apparent  increase  in  relevance  is  purely  artificial.   Hence  the  praxis 
of  considering  control  laws  to  be  only  stable  will  be  continued 
throughout  the  sequel. 

As  is  often  the  case,  in  this  theorem  practicality  is  the  price 
of  generality.   In  the  discussion  of  Theorem  2.3  it  was  observed  that 
the  conjugate  point  condition  was  not  very  satisfactory  as  a  sufficiency 
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condition  for  the  optimization.   The  same  remarks  apply  here:  no  con- 
ju[     points  will  exist  for  Q  positive  semidef inite;  then  the  conju- 
gate ;;oint  condition  need  only  be  examined  where  Q   is  not  positive 
semidef inite.   The  condition  that  Y  be  identical  for  equivalent  Q' s , 
part  (a) ,  can  also  be  simplified. 

The  next  lemma  provides  the  last  link  in  reformulating  the 
results  of  Theorem  3.3  into  a  workable  equivalence  relation  for 
performance  indices. 

Lemma  3.4 

For  an  arbitrary  real,  symmetric,  nth  order  matrix  Q, 
the  real  even  polynomial  resulting  from  the  quadratic  form 


Y(Q;u>)  =  S  (-jcu)QS(ju)) 


(15) 


where 


S(s)  = 


n-1 


can  be  determined  from  the  relation 

n  i-1 

Y(Q;u>)  =   E   [q.  .  +  2   E   (-l)Jq 
i=l   X1     j  =  l 


i-J ,i+J 


.]u> 


2(i-l) 


(16) 


(q.  .  =  0   for   i  <  1   or   j  >  n) . 


Proof 

Rewrite  the  vector  S  as  the  sum  of  two  orthogonal  vectors, 

0 


S(s)  =  a(s)  +  b(s)   a(s)  = 


b(s)  = 
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and  clearly 

S(-s)  =  a(s)  -  b(s) 

since  a  is  an  even  function  of  s  and  b  is  odd.   Substitution  of  this 

relation  for  S  into  (15)  leads  to 

Y(Q;s)  =  (a(s)  -  b(s))TQ(a(s)  +  b(s)) 

T       T        T       T 
=  a  Qa  -  b  Qb  +  (a  Qb  -  b  Qa)  . 

By  the  symmetry  of  Q  the  bilinear  terms  in  parentheses  total  zero. 

Since  a(s)  and  b(s)  have  odd  or  even  elements  respectively  zero, 

Y(Q;s)  contains  no  terms  with  q.  .  with  i  +  j  odd.   The  non-zero 

T 
off-diagonal  terras  resulting  from  a  Qa  have  i  and  j  both  even,  while 

T 
the  non-zero  off -diagonal  terms  from  b  Qb  have  i  and  j  both  odd. 

Hence,  Y  may  be  rewritten  as 

?(Q;i)  =   E  (-l)1+1q..s2<1-1>  +   S      q..si+j"2+  E       q.  .si+j"2 
1=1  i,j  odd  i,j  even 

which  reduces  easily  to  (16). 

This  lemma  is  important  in  its  own  right  in  that  it  provides 
an  algorithm  for  computation  of  Y  for  a  given  Q  without  the  necessity 
of  evaluating  the  quadratic  form.   Its  primary  value  however  is  that 
it  allows  the  equivalence  condition  of  Theorem  3.3  to  be  redefined  from 
the  invariance  of  a  polynomial  (Y)  to  specific  algebraic  constraints  on 
the  entries  of  the  matrices  in  question. 
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Corollary  3.3 

A  weighting  matrix  Q   is  equivalent  to  a  matrix  Q^   for 
the  system  of  Theorem  3.3  if  and  only  if  the  entries 
of  the  matrices  are  related  so  that  the  quantities 

p.  =  q. .  -  2q.    .    +  2q .  9  .  9  ~  •• •  < 17 > 

l    H      1-1,1+1      i-2,i+*2 

i  =  1,2,3.  .  .  ,n 

(q   =0  for  i  <  1  or  j  >  n) 
ij 

are  equal  for  q. .  taken  to  be  elements  Q  or  Q2  and 
parts  (b)  and  (9  of  Theorem  3.3  are  satisfied. 


Proof 

It  must  be  shown  that  the  relation  (17)  given  in  the  corollary 
is  equivalent  to  part  (a)  of  Theorem  3.3,  that  is,  if  it  is  satisfied, 
the  Y  polynomials  for  Q  and  Q  are  identical.   This  is  accomplished 
by  demonstrating  that  the  coefficients  of  the  respective  Y  polynomials 
are  coincident.   The  proof  follows  directly  from  determination  of  the 
coefficients  of  Y  from  (16)  which  are  then  related  to  the  pi  of  (17). 

With  additional  study  of  equations  (16)  and  (17)   a  somewhat 
startling  phenomenon  comes  to  light.   The  only  entries  of  a  given 
weighting  matrix  Q  which  influence  the  polynomial  Y(Q;w)  are  of  the  form 


qkk  and  qk-£,k+r 


which  excludes  any  element  q.  .  where  i  +  j  is  odd.   This  indicates 
that  approximately  half  of  the  elements  of  an  arbitrary  weighting  matrix 
are  irrelevant  (assuming  that  parts  b  and  c  of  Theorem  3.3  are  still 
satisfied)  with  respect  to  the  optimal  control  law.   Kalman  and  Englar 
[K3]  noted  from  the  structure  of  the  companion  matrix  canonical  form 
that  the  q.  .  elements  where  i+j  is  odd  are  "irrelevant"  without  identi- 
fying the  underlying  structural  relation  (17)  for  equivalent  Q's. 
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Their  suggestion  that  the  "  iri'elevant"  terms  in  a  given  weighting  matrix 
be  nulled  at  the  onset  of  the  optimization  procedure,  in  order  to  simplify 
computation,  is  potentially  hazardous.   There  is  the  possibility  that 
nulling  these  elements  will  alter  the  observability  qualities  of  the 
original  matrix  to  the  extent  that  an  unstable  plant  pole  is  obscured. 
Fortunately,  this  occurrence  appears  to  be  extremely  unlikely;  the 
author  was   unable  to  construct  a  matrix  which  behaved  in  this  fashion 
after  a  very  exhaustive  search.   It  seems  that  an  unobservable  matrix 
will  remain  so  after  the  i  +  j-odd  terms  have  been  struck  and  conversely 
an  observable  matrix  will  still  be  observable  after  these  terms  have 
been  removed.   A  very  convincing  heuristic  argument  can  be  made  in 
support  of  this  observation;  first,  however,  it  will  be  convenient  to 
state  two  lemmas. 

Lemma  3.5  (Gerschgorin,  1931  [G5]) 

All  the  eigenvalues  of  a  square  matrix  A  =  (a. .)  lie  in 
the  union  of  circular  regions, 

n 
|3ii  -  z|  £  2   |a   |,  i  =  1,2,3,. ...n, 

j  =  l    J 

of  the  complex  plane. 

This  is,  of  course,  the  touted  Gerschgorin  Circle  Theorem 
which  has  found  wide  application  in  the  numerical  eigenvalue  problem. 
A  laconic  proof  of  this  famous  theorem  can  be  found  in  Cull  en  [C2, 
p.  197].   Of  principal  immediate  interest,  however,  will  be  a  second 
lemma  which,  although  a  corollary  of  Gerschgorin' s  Theorem,  has  found 
prominence  as  a  separate  result. 
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Lemma  3. 6 

Diagonally  dominant  matrices  are  non-singular.   A  square 
matrix  A  is  said  to  be  diagonally  dominant  if 

n 

la   I  >  E   la   |,   for   i  =  1,2,3 n  (18) 

ii'    J=1   U 

where  a. .  is  the  (i,j)th  element  of  A;  that  is,  if 

the  diagonal  entries  are  larger  in  modulus  than  the  sum 

of  the  magnitudes  of  the  remaining  constituents  of  their 
respective  rows. 


The  justification  for  this  lemma  follows  easily  from  Lemma  3.5. 
If  a  matrix  is  diagonally  dominant,  the  region  of  permitted  eigenvalue 
locations  excludes  the  origin  and  the  matrix  is  consequently  non-singular. 

Consider  Figure  3.4  which  is  a  schematic  drawing  of  a  symmetric 
matrix  with  the  i  +  j  odd  terms  removed  and  the  X' s  representing  the 
remaining  entries. 

r* 

X   0    X    0    X  .  . 

Q  =    0    X    0    X    0  .  . 

X   0    X    0    X  .  . 


Figure  3.4   Sparse  Equivalent  Matrix 


It  is  clear  from  observation  of  the  figure  that  removing  the  indicated 
entries  tends  to  increase  the  dominance  of  the  diagonal  and  in  a  sense 
makes  it  "more  non-singular";  thus  it  would  be  reasonable  to  expect  an 
aggrandizement  of  observability  rather  than  a  deterioration. 

A  possible  second  concern  is  that  discarding  elements  in  the 
manner  described  may  destroy  the  positive  semidef initeness  of  the 
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weighting  matrix.   A  second  application  of  Gerschgorin* s  Theorem, 
this  time  in  the  form  of  Lemma  3.5,  reveals  that  the  likelihood  of 
degrading  the  positive  semidef initeness  of  a  Q  matrix  by  removing 
the  irrelevant  terms  is  extremely  small. 

If  the  matrix  is  positive  semidef inite,  the  union  of  circles 
which  form  the  permissible  regions  for  eigenvalues  must  include  parts 
of  the  right-half  complex  plane  (including,  perhaps,  the  origin);  and 
by  a  well-known  result  of  matrix  theory  [HI] ,  the  diagonal  elements  of 
a  positive  semidef inite  matrix,  hence  the  centers  of  these  circles, 
must  be  non-negative.   Then  as  the  off -diagonal  terms  are  removed,  as  in 
Figure  3.4,  it  is  apparent  from  relation  (18)  that  the  radii  of  the 
circles  which  may  contain  eigenvalues  are  reduced  and  the  eigenvalues 
will  consequently  be  restricted  to  fall  in  a  region  that  is,  if  anything, 
more  positive.   The  only  case  which  is  not  resolved  by  this  argument  is 
the  one  where  the  permissible  region  of  the  original  matrix  includes 
part  of  the  left  half-plane  and  reduction  of  the  radii  does  not  retrieve 
the  region  wholly  into  the  right  half-plane;  thus  admitting  the  possi- 
bility of  a  negative  eigenvalue  in  the  reduced  matrix.   In  any  case, 
the  migration  of  any  of  the  eigenvalues  of  a  Q  matrix  to  the  left,  as 
a  result  of  the  simplifying  operation,  appears  extremely  unlikely. 

In  general  bhe  observability  and  absence  of  conjugate  points 
for  the  sparser  matrix  must  be  tested.   However,  for  a  specific  but 
very  important  case,  both  conditions  (and  hence  equivalence)  can  be 
guaranteed  a  priori. 
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Theorem  3. 4 

If  a  weighting  matrix  Q  for  a  quadratic  performance 
index  (3)  operating  on  a  scalar  system  in  companion  matrix 
canonical  form  is  unity  rank,  i.e. , 

T 
hh   =  Qit 

where  h  is  a  vector,  and  the  resulting  optimal  control  law 

is  stable,  then  the  matrix  Q  resulting  from  nulling  the 

entries  q   of  Q  for  i+j  odd  is  equivalent  to  Q  . 


Proof 


Since  Q  is  unity  rank,  its  elements  can  easily  be  written  as 


functions  of  the  elements  of  the  factor  vector;  that  is, 

~  2 
ll 

2 

.  2 


Qx  =  hh  = 


hlh2 


hlh2 


V3 


h2h3 


hlh3 


h2h3 


hlh4 


2  4 


h3h4 


where  h.  is  the  ith  entry  of  h.   The  matrix  Q  constructed  from  Q  by 
discarding  the  elements  q.  .  where  i+j  is  odd  is 


Q2  = 


hlh3 


hlh3 


2  4 


It  is  easily  verified  that  a  rank  two  factor  of  Q,^  is 


H2  = 


\ 

0 

0 

h 

h3 

0 

Q2  =  H2H2 
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Part  (a)  of  Theorem  3.3,  i.e.,  YCQ-iO))  =  ¥(Q  ;cjo)  ,  is  satisfied 
by  virtue  of  Corollary  3.3  and  part  (c)  is  fulfilled  because  Q  is 
shown  to  be  positive  semidef inite  as  a  result  of  the  existence  of  the 
factor  H_.   Then  all  that  is  required  is  to  show  that  no  unstable 
poles  will  be  unobservable  by  H_  (part  b)  to  prove  Q  ~  Q  .   In  fact, 
an  even  stronger  condition  will  be  shown;  that  is,  the  only  plant  poles 
which  are  unobservable  through  Q   are  precisely  those  which  are  unobserv- 
able through  Q  .   This  will  effectively  satisfy  part  (b)  since  Q   is 
required  by  the  hypothesis  to  result  in  a  stable  control  law  and  hence 
cannot  fail  to  observe  any  unstable  plant  poles 

It  must  now  be  demonstrated  that  the  only  common  factors  of 
the  polynomial  entries  of  the  vector, 

H^SCs),  (19) 

are  also  factors  of  the  corresponding  polynomial  for  Q  , 

hTS(s)  . 

If  H_  is  thought  of  as  defining  two  outputs  to  the  system,  the  above 

requires  that  the  outputs  not  simultaneously  obscure  any  system  pole 

which  is  observable  through  Q  .   The  polynomials  formed  by  (19)  are 

2      4 
h   +  h  s   +  h  s   +  .  .  . 
1     o       5 

3    ,   5 
and  h  s  +  h  s   +  h  s   +  ...  . 

Z  4       b 

These  polynomials  can  be  recognized  as  the  result  of  "separating"  the 

T 
polynomial  h  S(s)  into  even  and  odd  functions  of  s.   This  is  a  common 

procedure  in    Neiwork  Theory  •;  in  which  the  polynomial  operators 

Ev  (•)  and  Od  (•)  are  often  defined  [T3]  to  simplify  notation. 
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In  this  notation 

F(s)  =  hTS(s)  =  Ev(F)  +  Od(F)  =  Ev(F)  +  sEv  (F) , 

where  Ev  (F)  is  Od(F)  with  an  s  factored  out.   Making  use  of  the  prop- 
o 

erty  of  even  polynomials  that  their  zeros  lie  symmetrically  about  the 

imaginary  axis  and  denoting  the  highest  even  power  coefficient  of  F 

as  c   and  the  highest  odd  power  coefficient  as  c  ,  F(s)  can  be 
e  o 


rewritten  as 


2  2 

F(s)  =  c   .TT  (s  +  a.)  +  s  c   .TT,  (s  +  b.). 
e  i=l        i       o  j=l        j 


In  the  above  ±  Va~~  and  ±  ^/b\    are  zeros  of  Ev(F)  and  Ev  (F>,  and 
1         J  o 

the  upper  limits  on  the  products  must  be  chosen  with  regard  to  whether 

the  dimension  of  the  system  is  even  or  odd.   Now  if  Ev(F)  and  Ev  (F) 

possess  a  common  factor,  then  a,  =  b  for  some  k  and  m  and 

k    m 

F(s)  =  (s2  +  a,  )  c   TT   (s2+a.)  +  c   TT  (s2+b  )] 
k  U  e  i=l      !     o  j  J 

which  is  clearly  a  factor  of  F(s)  as  well,  and  the  theorem  is  proven. 

This  proof  of  the  observability  of  the  rank  two  Q  will  be  use- 
ful in  the  next  theorem.   A  simpler  analysis  follows  easily  from 
consideration  of  the  composition  of  the  two  outputs  defined  by  H  . 
Summing  these  outputs  results  in  a  single  output  identical  to  the  one 
defined  by  h,  thus  any  response  which  is  observable  from  h  is  also 
observable  from  H  . 

This  theorem  guarantees  that  if  a  system  and  an  optimal  control 
minimize  a  single  output  in  the  mean-square  sense  (h  is  unity  rank)  , 
then  with  little  effort  a  pair  of  outputs  can  be  defined  which  are  also 
minimized.   The  real  value  of  this  result  is,  however,  that  it  allows 
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enormous  simplification  of  the  unity  rank  performance  index  at  the  onset 
of  the  problem,  with  attendant  savings  in  numerical  quality  and  quantity. 
This  result  will  also  find  application  in  the  next  chapter  in  another 
context. 

A  special  case  of  Theorem  3.4  reveals  an  interesting  structure 
which  is  well  worth  recording. 


Corollary  3.4 

If  matrix  Q..  of  Theorem  3.4  is  such  that  its  factor  h 
forms  a  polynomial  hTS(s)  which  has  all  zeros  with  non- 
positive  real  parts;  then  the  factor  H2  of  matrix  Q2  con- 
structed as  described  in  the  proof  will  form  two  polynomials 


F1(s) 
F2(s) 


H^SCs) 


having  simple  zeros  (except  possibly  at  the  origin)  which 
are  restricted  to  lie  on  the  imaginary  axis  where  they 
occur  in  conjugate  pairs  and  alternate  with  each  other. 


Proof 

The  proof  only  requires  that  the  root  location  property 
described  be  demonstrated.   The  polynomials  of  interest  are  easily 
recognized  as 

F1(s)  =  Ev(F) 
F2(s)  =  Od(F) 
and    F(s)  =  hTS(s). 

It  is  a  well-known  theorem  of  Network  Theory  that  the 
Ev(F)   and  Od(F)   functions  of  a  Hurwitz  polynomial,  F,  have  the  root 
partitioning  property  recounted  in  the  corollary  [G6],   This  phenom- 
enon is  referred  to  as  the  Alternation  [G6]  or  Separation  TT3]  property. 
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A  complex  plane  diagram  of  the  respective  zero  locations  of  a  typical 
case  is  useful  to  visualize  this  result. 


■*■ 


Re-» 
o 
i 
x 


o  -  zeros  of  Ev(F) 
x  -  zeros  of  Od(F) 


i  F  -  a  Hurwitz  polynomial 


Figure  3.5   Zero  Locations  of  Ev(F)  and  Od(F) 


This  corollary  implies  more  than  appears  in  the  hypothesis. 

T 
By  Theorem  3.2  a  unity  rank  Q,  with  h  S(s)  a  Hurwitz  polynomial,  can 

always  be  constructed  if  Y(u))  is  non-negative,  then  the  sparse  matrix 

of  the  corollary  can  always  be  constructed  as  well. 


3. 6   Resume  of  Y-Invariant  Matrices 

Because  of  the  necessity  of  investigating  observability  and 
the  testing  for  the  absence  of  conjugate  points  much  of  the  heuristic 
appeal  of  the  form  of  matrices  invariant  under  the  operator  Y(Q;u>) 
is  obscured.   In  order  to  distinguish  between  matrices  equivalent  in 
the  sense  of  Theorem  3.3  and  matrices  which  result  in  identical  Y' s  , 
without  resorting  to  cumbersome  phraseology,  the  latter  will  be 
referred  to  as  "Y-invariant  matrices. "  The  conjugate  point  condition 
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is  satisfied  a  priori  for  a  very  large  class  of  weighting  matrices 
(i.e. ,  positive  semidef inite)  and  the  observability  restriction  has 
been  shown  to  be  a  relatively  innocuous  constraint  between  ^'-invariant 
matrices.   Then  a  great  deal  is  to  be  gained  from  a  review  of  the 
structure  of  invariant  matrices. 

General  Structure 


Corollary  3.3  reveals  that  two  matrices  are  Y-invariant  if  and 
only  if  their  elements  are  such  that  the  n-tuple 


Pi  =  qii  "  2qi-l,i+l+  2qi-2,i+2 


(q  .  =  0   for   i  <  1   or   i  >  n) 


i  =  1,2,  .  .  .  ,n 


is  equal  for  both  matrices.   For  instance,  in  the  third  order  case  this 
can  be  interpreted  as  meaning  that  Y  will  remain  invariant  if  a  quantity 


is  added  to  q     and  q     and  twice  that  quantity  is  added  to  q 


22 


In  general ,  matrices  which  are  Y-invariants  of  a  given  matrix  may  be  con- 
structed by  manipulating  the  elements  on  "diagonals"  running  from 
lower  left  to  upper  right.   This  can  be  represented  pictorially  as  in 
Figure  3.6. 


n-1 


n-1 


n-1 


Figure  3.6   Structure  of  Y-Invariant  Matrices 
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The  integer  entries  of  the  illustration  indicate  which  terra  of  the 
n-tuple  is  affected  by  the  element  of  a  Q  matrix  which  would  be  in  that 
location;  the  x1 s  are,  of  course,  entries  which  influence  no  term. 
From  this  it  is  easy  to  see  that  matrices  Y-invariant  <to  a  given  matrix 
can  be  identified  by  inspection. 

Spectral  Factorization 

A  special  case  which  is  of  considerable  interest  is  the  Q 
formed  from  the  coefficients  of  the  Hurwitz  spectral  factor  of  Y(uu)  . 
It  is  a  form  which  can  always  be  constructed  when  Y(cu)  is  non-negative 
and  always  meets  the  observability  and  conjugate  point  criteria.   When 
the  i  +  j  odd  terms  are  zeroed  the  resulting  matrix  is  always  equivalent 
and  possesses  the  root  partitioning  property  of  Corollary  3.4.   This 
form  is  theoretically  straightforward  to  compute  but  in  practice  is 
among  the  most  difficult. 

Diagonal 

Another  special  case  which  is  appealing  in  its  simplicity  is 

the  diagonal  Q.   This  matrix  is  formed  by  placing  the  coefficients  of 

2 
Y(uu  )  along  the  diagonal  with  the  constant  term  first.  A  Y-invariant 

diagonal  matrix  also  can  always  be  formed;  no  observability  difficulties 

will  be  encountered  if  the  terras  are  all  non-zero  and  the  conjugate 

point  condition  will  be  satisfied  if  the  terms  are  non-negative. 

A  simple  example  helps  to  demonstrate  how  this  all  fits  together. 


GO 


Example  3. 3 

Consider  the  real,  even  polynomial, 


4      2 
T(ud)  =  uj  +  2<D  +1, 


and  some  members  of  the  set  of  ( Y-invariant)  matrices, 

{q|y(Q;cu)  =  r(uo)}  . 
i)   Diagonal 

The  diagonal  Y-invariant  matrix  is  simply, 


Ql  = 


1  0  0 
0  2  0 
0     0      1 


which  is  non-singular  (thus  observable)  and  also  positive  definite. 
ii)   Unity  Rank 

The  roots  of  r(ai)  are  sketched  in  the  complex  plane^'diagram 
below. 


-1 


4  1 


If    the   left   half-plane   zeros    are   chosen,    the  resulting   unity   rank 


T 
=    I£H      is 


*2  = 


12  1 
2  4  2 
12         1 


[1         2         1] 


and  zeroing  the  i+j  odd  terms  results  in 
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-Q3  = 


~ion 

0         4        0 
10         1 


'i    on 

0  2 

1  0 


n     o     n 

0         2         0 


This  matrix  is   clearly   rank   two   and   can   easily   be  shown  a ¥- invariant 
of   Q     by    adding    -1   to   q, _    and  q        and   twice   that    (-2)    to   q      .      The 

1  Id  Jl  ^^ 

T 
zeros  of  H  S(s)  for  Q  are,  respectively,  ±j  for  the  first  output  and  0 

for  the  second,  demonstrating  the  separation  property  of  Corollary  3.4. 

If  the  factor  is  composed  of  one  zero  from  each  half-plane,  the 

resulting  Q  is 


Q  = 
*4 


0    -1 


0    0 


which  is  again  clearly  an  invariant  of  Q  . 

Note  that  in  every  Q  matrix  gener   ed  in  this  example  the  q 


11 


and  q    entries  remain  fixed.   From  Corollary  3.3  it  is  obvious  that 
nn 

this  must  always  be  the  case  and  that  these  terms  are,  respectively, 
the  first  (constant  term)  and  last  coefficients  of  Y(uu) .   The  constancy 
of  these  elements  has  a  very  interesting  physical  interpretation  which 
will  be  discussed  in  the  next  chapter. 


3.7   Summary 

The  introduction  to  this  chapter  alluded  to  a  concern  which  has 
been  often  expressed,  with  varying  degrees  of  vehemence,  by  those  well- 
versed  in  classical  design  techniques;  that  is,  the  minimality  of  some 
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performance  measure  is  seldom  directly  relevant  to  a  practical  design 
problem.   For  the  quadratic  performance  index  and  a  linear  system, 
however,  it  was  shown  that  resulting  optimal  system  configurations 
possess  characteristics  which  are  very  much  in  the  spirit  of  good 
system  design.   The  minimum  optimal  system  stability  margins  of 
60°  in  phase  and  6  db  in  gain  compare  favorably  with  those  encountered 
in  practice. 

The  problem  of  actually  constructing  a  performance  index  which 
is  minimized  by  an  optimal  system  was  found  to  lead,  with  due  consider- 
ation for  certain  cancellation  and  existence  difficulties,  to  an  equiv- 
alence class  of  weighting  matrices  which  enjoy  some  extremely  enlight- 
ening and  useful  mutual  properties.   When  a  single  member  of  an  equiva- 
lence class  has  been  determined. the  remainder  of  the  class  can  be 
constructed  with  relative  ease. 


CHAPTER  IV 
THE  INVERSE  PROBLEM  AND  LINEAR  REGULATOR  DESIGN 

4.1   Introduction 

The  previous  chapters  have  dealt  with  optimal  systems ,  their 
properties  and  construction.   It  is  essential  that  these  concepts  be 
well  in  hand  if  the  techniques  of  optimal  control  theory  are  to  be 
applied  to  design  problems  where  a  relevant  performance  measure  is  not 
readily  discernible.   In  general,  this  is  the  case  when  a  linear  regu- 
lator is  to  be  designed  using  optimal  control  theory.   This  chapter 
will  investigate  how  optimal  control  theory  and  classical  techniques 
may  be  used  to  complement  one  another  in  practical  design  problems. 

For  the  purposes  of  this  chapter  it  is  helpful  to  think  of 
classical  synthesis  techniques  for  linear  control  systems  in  terms  of 
Figure  4.1.   The  problem  originates  with  a  linear  system  and  a  number 
of  performance  specifications  that  the  compensated  system  is  to  satisfy 
(node  A).   The  desired  result  is  a  realization  (node  C)  which  meets  the 
performance  criteria  and  contains  a  compensator  which  is  satisfactory 
from  the  standpoint  of  practicality  constraints,  such  as  realizability 
and  noisy  or  incomplete  measurements. 

One  way  of  approaching  the  problem  is  to  prescribe  pole  loca- 
tions which  will  meet  the  performance  specifications  (node  B)  and  then 
to  construct  a  compensator  which  will  place  the  plant  poles  in  approx- 
imately the  desired  locations  while  not  contributing  significantly 
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to  system  response  (p;ith  II).   This  is  a  quite  prevalent  philosophy, 
although  it  is  often  obscured  by  the  specific  design  procedure  used. 
Path  I  can  be  thought  of  as  classical  synthesis  and  analysis  procedures 
used  iteratively  to  arrive  at  these  pole  locations. 

Specifications 
A 


Realization 


Pole 
Locations 


Classical  Design  Procedure 
I.   Prescribe  pole  locations 
II.   Design  compensator 
III.   Direct  design 


Alternatives 

1.  Compute  performance  index 

a.  from  specifications 

b.  from  pole  locations 

2.  Optimal  control  problem 

3.  Automated  optimal  compensator 
design 


Figure  4. 1   Control  System  Design  Procedures 


A  second  approach  is  to  determine  at  the  onset  the  form  of  the 
compensator  required  and  manipulate  its  parameters  to  arrive  directly 
at  a  realization  (path  III).   Well-known  variations  of  this  approach 
are  Evans'  root  locus  and  lead-lag  design  [E2] .   These  techniques  have 
the  decided  advantage  of  defining  a  clear-cut  way  to  proceed  but  it  may 
be  difficult  to  accommodate  noisy  measurements. 

The  procedures  of  this  chapter  provide  alternative  paths  in  the 
design  scheme  as  illustrated  in  Figure  4. 1  by  the  Arabic  numbered 
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branches.   This  provides  for  much  greater  flexibility  in  how  the  problem 
is  attacked  and  employs  the  power  of  optimal  control  theory  to  relieve 
some  of  the  procedural  or  computational  burden  in  some  or  all  phases 
of  the  design. 

In  some  cases  it  will  be  possible  to  define  a  performance  index 
directly  from  classical  specifications  (path  la)  which  will  be  minimized 
by  a  control  law  which  satisfies  the  specifications.   The  specific  con- 
trol law  can  then  easily  be  computed  (path  2)  or  the  corresponding 
compensators  can  be  designed  (path  3) ,  using  Kalman-Bucy  filter  theory 
|"K8,K9]  or  one  of  the  new  techniques  for  automated  compensator  design 
[P1.P2]. 

As  indicated  by  path  lb  of  the  figure,  it  is  possible  to  compute 
a  performance  index  corresponding  to  a  closed-loop  (optimal)  pole  con- 
figuration and  the  remainder  of  the  synthesis  can  be  completed  through 
the  use  of  compensator  design  techniques  discussed  in  the  last  paragraph. 

Figure  4.1  does  not  reveal  some  of  the  additional  flexibility 
allowed  by  these  procedures.   For  instance,  a  performance  index  may 
be  specified  which  meets  only  a  subset  of  the  specifications  and  may 
then  be  used  as  a  basis  for  design  iteration.   The  techniques  also 
provide  for  the  design  of  sampled-data  compensators  (Section  4.4). 

4. 2   Pole  Placement  by  Performance 
Index  Designation 

If  a  set  of  closed-loop  system  poles  are  proposed  as  a  prelim- 
inary or  final  design  configuration,  the  prescription  of  a  performance 
index  which  is  minimized  by  this  configuration  (path  lb  of  Figure  4.1) 
may  be  of  considerable  value  in  the  completion  or  refinement  of  the 
design.   Figure  4.2  illustrates  such  a  case. 
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n — n  —  x — x- 


x  -  plant  poles 

□  -  proposed  closed-loop  poles 


Figure  4.2   Proposed  Closed-loop  Pole  Configuration 


This  configuration  may  have  resulted  from  a  root-locus  type 
analysis,  the  exact  scheme  used  is  not  germane  to  the  present  discus- 
sion; the  important  point  is  that  it  is  not  bound  to  any  technique  or 
form  of  compensation.   The  construction  of  a  performance  index  which 

is  minimized  by  this  configuration  (if  possible)  provides  at  most 

2  1 

approximately  n  /4  nominal  parameters   which  may  be  varied  to  perturb 

the  design  or  used  to  design  the  required  compensator  (path  3  of 
Figure  4. 1) . 

The  problem  of  computing  a  quadratic  performance  index  which 
is  minimized  by  specified  closed-loop  pole  locations  for  a  given 
plant  is  basically  the  problem  entertained  in  great  detail  in  the 
preceding  chapter.   This  section  seeks  to  deal  with  this  problem  on 
a  more  pragmatic  basis. 

The  plant  to  be  considered  is  a  completely  controllable, 
constant,  linear  system  taken  without  loss  of  generality  to  be  in 
companion  matrix  canonical  form  and  the  performance  index  is  the 
infinite  final  time  case  of  the  last  chapter. 


The  precise  number  of  pertinent  parameters  will  be  given  later. 
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The  stable  closed-loop  pole  configuration  of  Figure  4.2  will  be 
optimal  with  a  positive  semidefinite  weighting  matrix  Q  if  and  only  if 

cp(io)  =  |<p  (ju)|2  -  |cp(ju>)|   ^  0    for  all  real  cu  (1) 

where  cp(s)  and  9  (s)  are  the  normalized  (highest  coefficient  is  unity) 
open-  and  closed-loop  characteristic  polynomials.   Although  (1)  is  a 
succinct  optimal ity  criterion,  it  is  not  obvious  how  best  to  proceed 
to  test  a  given  polynomial  (1)  for  non-negativity. 

The  necessity  for  testing  the  sign  semidef initeness  of  real 
even  polynomials  arises  in  many  other  applications,  for  instance,  tests 

for  positive  reality  in  Network  Theory  [VI].   It  can  be  shown  by 

2 

factoring  Y(cu  )  or  through  the  use  of  Sturm's  Theorem  that  (1)  is 

2 
equivalent  to  requiring  that  Y(uu  )  have  no  positive  real  roots  of  odd 

multiplicity  [VI,  p.  106]. 

There  is  still  not  a  computationally  satisfactory  approach 

2 

evident.   Calculation  of  the  zeros  of  Y(u>  ),  if  the  order  is  large, 

is  difficult  numerically  and  best  avoided  as  a  test.   Recent  results 
by  Siljak  [S3]  and  Karmarkar  [K10]  extend  the  interpretation  of  the 
sign  changes  in  the  first  column  of  the  Routh  table  [G4]  to  test  for 
positive  real  zeros  of  odd  multiplicity.   This  test  is  probably  the 

only  alternative  currently  available  to  the  computation  of  the  zeros 

2 

of  Y(cu  )  as  a  necessary  and  sufficient  test  for  non-negativity. 

Appendix  C  outlines  a  numerical  implementation  (in  Fortran  IV) 
of  a  modification  of  Siljak' s  method,  which  is  both  accurate  and  effi- 
cient; it  is  believed  to  be  the  only  such  program  in  existence. 
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It  should  be  reiterated  that  condition  (1)  applies  only  to 
optimality  with  respect  to  a  positive  semidefinite  Q;  a  system  with 
Y(uj)  }£   0  may  well  be  optimal  for  a  sign  indefinite  Q.   In  such  a  case 
the  properties  conditioned  on  Q  being  positive  semidefinite  (Section 
3.3)  would  not  in  general  hold  but  the  Q  equivalence  relations  and 
remarks  concerning  the  observability  requirements  for  Q  (Chapter  III)do. 

If  criterion  (1)  is  satisfied,  a  positive  semidefinite  weighting 
matrix  can  be  generated,  in  fact,  an  entire  equivalence  class  of  posi- 
tive semidefinite  Q's.   In  general  the  initial  member  of  the  equivalence 
class  must  be  computed  by  spectral  factorization.   This  operation,  like 
the  test  for  non-negativity,  is  not  a  numerically  trivial  one.   The 
spectral  factorization  of  Y(tw)  ,  if  approached  naively,  consists  of 
determining  the  2(n-l)  roots  of  Y,  separating  them  by  the  sign  of  their 
real  parts,  constructing  the  factor  composed  of  the  zeros  with  negative 
real  parts  and  multiplying  it  by  the  appropriate  constant.   This  pro- 
cedure is  entirely  unsatisfactory.   It  is  well  known  that  the  confidence 
in  approximation  for  root  locations  generally  decreases  drastically  for 
polynomials  of  large  order;  this  coupled  with  the  error  induced  by 
constructing  the  spectral  factor  from  its  roots  makes  this  procedure 
numerically  hazardous. 

An  obvious  alternative  is  to  reduce  the  order  of  Y(u))  by  sub- 

2 

stituting   c  =  id   and  computing  the  roots  of  an  n  -  1  order  polynomial 

which  are  the  squares  of  the  actual  roots  of  interest.   The  roots  of 
Y(cu)  with  negative  real  parts  can  be  obtained  directly.   A  refinement 
of  this  procedure  is  to  compute  quadratic  factors  of  Y(a)  only  and 
from  these  factors,  through  some  rather  intricate  logic,  garner  the 
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related  left  half-plane  quadratic  factors  of  Y(co) .   This  serves  to 
reduce  the  total  number  of  numerical  manipulations,  saving  computing 
time  and  decreasing  the  sources  of  error  propagation.   The  implementa- 
tion of  the  spectral  factorization  algorithm  discussed  in  the  appendix 
takes  this  approach,  utilizing  an  efficient  technique  due  to  Bairstow 
[Kll],  for  the  approximation  of  quadratic  factors.   Bairstow' s  iter- 
ation has  been  shown  to  have  a  rapid  rate  of  convergence,  although  it 
is  somewhat  more  sensitive  to  starting  values  than  competing  methods 

[Kll,  p.  lOlff]. 

One  objection  to  this  procedure  is  that  the  roots  of  Y(uj)  are 
not  directly  available  as  a  secondary  test  for  non-negativity.   The 
test  for  non-negativity  described  earlier  seems  wholly  adequate  and 
any  sacrifice  of  computational  efficiency  in  the  spectral  factorization 
in  favor  of  a  redundant  test  does  not  appear  justifiable. 

The  flow  diagram  of  Figure  4.3  summarizes  concisely  the  steps 
to  be  taken  in  the  generation  of  a  quadratic  performance  index  which 
is  minimized  by  a  given  plant  and  stable  proposed  (optimal)  closed-loop 
pole  configuration.   The  first  step  (T)  is  to  compute  the  closed-loop 
and  plant  normalized  characteristic  polynomials  (if  not  already  known). 
The  magnitude-square  of  the  open-  and  closed-loop  characteristic  poly- 
nomials are  next  computed;  this  is  another  operation  where  the  obvious 
procedure,  i.e.,  multiplication  of  the  respective  polynomials  by  their 
complex  conjugates,  is  not  the  best  choice.   Since  the  magnitude-square 
polynomials  are  even,  half  of  the  coefficients  are  zero.   The  technique 
outlined  in  Appendix  C  makes  use  of  this  structure  by  computing  only 
the  non-zero  terms  and  in  a  way  which  obviates  the  use  of  complex 
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Figure  4.3   Procedure  for  Specification  of  a  Performance 
Index  from  Proposed  Pole  Locations 


Numbers  in  circles  are  referred  to  in  the  text. 
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arithmetic.   When  the  difference  of  the  magnitude-square  polynomials 
is  taken  (2),    the  resulting  Y(uu)  is  an  even  polynomial  of  order  2(n-l) 
due  to  the  normalization  of  the  highest  power  coefficients  of  ep  and  9 
to  one. 

A  decision  based  on  the  non-negativity  of  Y  is  made  in  block(3\ 
If  Y  is  non-negative  the  system  is  clearly  optimal  and  the  specifica- 
tion of  the  Q1 s  may  proceed.   The  unity-rank  Q  is  computed  from  the 
left  half-plane  factor  of  Y  Q*),    using  the  techniques  for  spectral 
factorization  discussed  earlier,  which  further  insures  the  observabil- 
ity conditions   are  satisfied,  and  the  equivalence  class  of  positive 
semidefinite  Q' s  is  immediately  availablef5\ 

If  Y  is  not  non-negative  the  procedure  is  no  longer  straight- 
forward and  positive  results  are  not  guaranteed.   An  initial  Q  is 
constructed  so  that 

Y(Q  ;cu)  =  Y(uj)  , 

and  all  the  unstable  plant  poles  are  observable/ 6j.      A  reasonable 

2 
choice  for  an  initial  Q  is  the  diagonal  case  (coefficients  of  Y(uu  ) 

on  diagonal) ;  this  is  the  easiest  to  generate  and  a  desirable  form 

for  the  answer.   This  choice  of  Q  is  inserted  in  the  Riccati  equation 

for  the  system  in  companion  matrix  canonical  form . 

P(t)  =  -FTP(t)  -P(t)F  +  P(t)ggTP(t)  -Q,  P(t)=0.  (2) 

The  Riccati  equation  is  then  integrated  numerically  until 

||  P(t)||  £  e,  (3) 

where  e  is  some  prescribed  non-negative  bound  or  until  a  conjugate 
point  intervenes  (the  solution  becomes  unbounded)(7 j.      The  non-linear 
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nature  of  the  Riccati  equation,  which  makes  this  sort  of  test  necessary, 
also  helps  insure  the  success  of  the  test,  since  it  may  be  subject  to 
the  phenomenon  of  finite  escape  time  [K2].   That  is,  as  opposed  to  the 
solutions  to  linear  differential  equations  which  may  be  unbounded  only 
in  the  limit,  the  Riccati  equation  may  become  unbounded  in  a  finite  time. 
Experience  has  shown  that  when  a  conjugate  point  does  exist  (2)  generally 
becomes  unbounded  quite  rapidly  (in  terms  of  the  number  of  integration 
steps)  and  conversely  when  a  conjugate  point  is  not  present  steady- 
state  (3)  is  reached  promptly. 

A  third  possibility  is  that  (2)  has  neither  a  conjugate  point 
nor  a  steady-state  solution;  that  is,  the  solution  is  oscillatory. 
It  is  clear  that  if  (2)  is  oscillatory  it  must  also  be  periodic  and 
hence  the  cost  and  the  system  response  are  periodic.   This  would 
contradict  the  supposition  that  the  proposed  system  design  was  time- 
invariant  and  stable.   Then,  if  a  periodic  solution  to  the  Riccati 
equation  occurs,  it  must  be  due  to  numerical  errors  and  discarded  in  the 
same  manner  as  a  solution  with  a  conjugate  point. 

A  Q.  which  results  in  a  solution  of  (2)  which  diverges  obviously 

must  be  discarded  but  it  does  not,  in  itself,  indicate  that  the  system 

is  not  optimal.   The  procedure  branches  to  block  {8 J   if  a  conjugate 

point  is  present  and  Q    which  is  Y-invariant  to  Q.  is  computed  and 
J  +  l  J 


th 


e  test  for  conjugate  points  [7j   is  repeated.   It  is  difficult  to 
conceive  of  an  algorithmic  scheme  for  "improving"  Q.  in  block  (  8  \ 
In  fact,  if  this  technique  were  known,  one  could  apply  it  iteratively 
until  the  "best"  Q  is  obtained  and  the  subsequent  test  for  conjugate 
points  would  be  a  general  necessary  and  sufficient  test  for  the 
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optimality  of  a  control  law   An  interactive  strategy  is  ideal,  ¥- 
invariant  Q' s  can  be  generated  which  not  only  are  more  likely  to  con- 
verge but  also  possess  a  desirable  structure  for  the  specific  applica- 
tion. 

The  algorithm  to  be  described  seeks  to  generate  anY-in^ariant 

Q    from  a  Q  which  fails  the  conjugate  point  condition  (7)  by  shift- 
j+1         J  v-y 

ing  the  negative  eigenvalues  of  Q.  generally  to  the  right  while  main- 
taining as  much  simplicity  as  possible  in  Q.   .   This  is  done  by  test- 

2 
ing  the  diagonal  terms  q. .  of  Q.  for   1  <  i  <  n  progressively   until 

a  negative  entry  is  detected  which  is  then  deleted,  using  equation  (17) 

of  Chapter  III. 

Suppose  that  the  minor  shown  below  is  a  3  x  3  principal  minor, 

containing  the  negative  diagonal  entry  b. 


Because  the  first  Q  (Q.)  was  diagonal  and  since  the  method  proceeds 
progressively  down  the  diagonal,  the  remainder  of  the  rows  shown  are 
null  except  the  first  which  may  have  an  off -diagonal  entry  -d/2  placed 
by  the  preceding  iteration.   By  the  results  of  the  last  chapter, 
altering  this  minor  as  shown  below: 


As  noted  in  the  preceding  chapter,  the  first  and  last  terms  on 
the  diagonal  cannot  be  altered  and  invariance  maintained. 
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leaves  the  V  of  the  resulting  matrix  invariant.   Each  time  a  Q  •  fails  the 
conjugate  point  test  the  described  procedure  removes  an  additional 
negative  diagonal  entry  until  all  have  been  eliminated.   In  this 

manner  the  scheme  allows  k+1  iterations  where  k  is  the  number  of  nega- 

2 

tive  coefficients  of  Y(uu  )  (or  the  number  of  negative  diagonal  terms) 

before  the  system  is  discarded.   The  procedure  is  executed  iteratively, 
rather  than  eliminating  all  of  the  negative  terms  in  one  step,  with 
the  hope  of  determining  as  simple  (diagonal)  a  Q  as  possible  for  which 
the  system  is  optimal. 

At  this  point  the  justification  for  such  a  procedure  is 
probably  not  clear   The  object  of  the  algorithm  is  to  generate 
a  sequence  of  invariant  Q' s  whose  positive  diagonal  elements  dominate 
their  respective  rows.   This  will  hopefully  result  in  a  shifting  of 
the  eigenvalues  to  the  right.   Central  to  this  argument  is  the 
Gerschgorin  Circle  Theorem  which  was  recorded  as  Lemma  3.5  in  the 
previous  chapter. 

Since  the  rows  of  any  Q.  generated  by  the  described  technique 
have  at  most  two  off-diagonal  elements,  each  being  half  of  a  negative 
diagonal  term,  the  eigenvalues  of  the  modified  matrices  will  be 
increasingly  restricted  to  be  "close"  to  the  right  half-plane  if  the 
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3 
diagonal  elements  are  approximately  of  the  same  magnitude.    There  is 

no  hope  that  the  matrices  will  eventually  become  positive  semidefinite 

with  repeated  iterations  (this  would  require  Y(uu)  ^  0  by  Theorem  3.2); 

however,  it  is  reasonable  to  expect  the  likelihood  of  passing  the 

conjugate  point  condition  to  increase  if  the  eigenvalues  are  shifted 

to  the  right.   An  example  will  help  to  clarify  the  details. 

Example  4. 1 

Consider  the  sequence  of  Q1 s  which  would  be  generated  by  this 

method  for 

ip     in     ft     fi     4     2 
Y(cu)  -7u)        -6u)   -5a)+3u)+2tu-4u)+l. 

This  Y  is  clearly  not  non-negative,  e.g. ,  Y(l)  =  -  2;  then  there  is  no 
possibility  of  constructing  a  positive  semidefinite  Q  for  it.   The 
procedure  beginning  with  block  C6J   of  the  flow  chart  of  Figure  4.3 
must  be  called  upon.   Since  there  are  three  negative  coefficients 
of  Y  there  will  be  at  most. four  ^-invariant  Q1  s  generated;  Figure  4.4 

illustrates  these  iterations.   The  first,  Q  ,  obviously  has  eigenvalues 

2 
which  are  simply  the  coefficients  of  t(cu  )  .   The  first  iteration  removes 

the  first  negative  entry  on  the  diagonal  (-4)  by  "splitting"  it  into 
the  two  off -diagonal  terms  resulting  in  Q^.   Four  of  the  eigenvalues 
remain  unchanged,  the  one  that  was  at  -4  is  moved  to  the  origin  and 
the  remaining  two  have  been  "smeared"  by  the  interpretation  given  by 
Gerschgorin' s  Theorem.   The  eigenvalue  which  was  at  1  may  now  be  any- 
where within  the  interval  (-1.+3)  and  the  eigenvalue  at  2  may  now  lie 
in  the  region  (0,+4)  as  shown  in  the  plot. 

3 

Actually  the  right  half-line,  since  the  matrices  under  con- 
sideration are  symmetric  and  consequently  have  only  real  eigenvalues. 


12      10      8      6      4      2 

Y(u>)  =7o)   -  6a>   -5cu+3a)+2(ju-4a)+l 
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With  subsequent  iterations  the  original  eigenvalues  are  spread 
into  permitted  regions  which  are  generally  to  the  right.   Due  to  the 
structure  of  such  matrices  the  actual  eigenvalues  are  easily  computed 
and  included  in  the  plot.   The  process  lives  up  to  expectation,  the 
negative  eigenvalues  are  moved  distinctly  to  the  right  and  the  likeli- 
hood of  passing  the  conjugate  point  condition  (esp.  for  Q  )  is  con- 
siderably enhanced. 

Note  that  the  rank  of  Q  for  j  >  1  is  reduced  by  at  least  one. 
This  causes  concern  that  some  observability  may  be  lost  in  the  process, 
with  the  result  that  an  unstable  plant  pole  becomes  unobservable  (a 
stable  plant  pole  becoming  unobservable  is,  of  course,  inconsequential). 
Although  this  may  indeed  be  the  case,  it  will  be  shown  later  that  suf- 
ficient safeguards  exist  to  prevent  it  from  passing  unnoticed. 

Within  the  framework  of  the  process  outlined  there  exist  many 
variations  which  may  be  of  value  in  certain  applications.   For  instance, 
it  may  be  advantageous  to  dispatch  the  negative  diagonal  terms  of 
largest  modulus  first  rathe:*'  than  approaching  them  progressively  along 
the  diagonal.   This  modification  applied  to  Example  4.1  leads  to  a 
somewhat  greater  "improvement"  in  the  second  iteration  than  the  method 
originally  employed,  but  the  results  of  the  final  iteration  are,  of 
course,  identical.   The  important  observation  to  make  is  that  no  single 
variation  will  be  best  in  every  case  and  that  employment  of  good  judg- 
ment at  this  point  in  the  scheme,  rather  than  following  a  strictly 
algorithmic  approach,  will  probably  be  rewarded  with  simpler  resulting 
Q's.   It  is  for  this  reason  that  the  suggestion  was  made  earlier  to 
utilize  an  interactive  strategy  in  a  computer  implementation  of  this 
procedure. 
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If  the  iterative  technique  of  modifying  Q.  and  testing  for 
conjugate  points  succeeds  in  obtaining  a  Q  for  which  the  Riccati  equa- 
tion does  not  diverge,  the  next  step  is  to  determine  if  this  Q  is 
actually  a  solution  to  the  inverse  problem.   This  is  necessary  to  guard 
against  gross  numerical  errors  in  the  integration  of  the  Riccati 
equation  falsely  indicating  the  absence  of  conjugate  points,  and,  as 
mentioned  earlier,  to  insure  that  a  detrimental  loss  of  observability 
from  Q  has  not  occurred. 

The  steady-state  solution,  P,  to  the  Riccati  equation  (2)  is 
available  from  the  conjugate  point  test  (block  (g\    and  provides  this 
check.   The  feedback  control  law, 

k  =  Pg, 
is  simply  the  last  column  (and  row)  of  P.   By  the  structure  of  the 
companion  matrix  canonical  form  if 

cp(s)  =  s  +  a  s    +  .  .  .  +  a0s  +  a_  , 
n  2     1' 

then    cp  (s)  =  s   +  (a  +k  )s    +  ...  +  (a„+k0)s  +  (a  +k  ), 
K  n   n  2   2       11 

where  k.  is  the  ith  entry  of  k.   The  coefficients  of  cp   computed  by 

substituting  p.    for  k.  can  be  compared  to  the  correct  (original)  values. 

This  test  is  represented  by  block  no)    of  Figure  4.3. 

In  passing,  it  should  be  noted  that  the  process  for  generating 

successive  Q' s  outlined  here  fails  to  make  full  use  of  the  relation 

of  the  elements  of  invariant  matrices  (Corollary  3.3)  in  that  only  the 

terms  on  the  diagonal  and  second  superior  diagonal  are  affected  (i.e., 

q.  .  ,  q.  ,,  .  _  and  q.  „  .  „  when  q  .  <  0)  .   This  is  because  the  process 
n    i+2,i-2       i-2,i+2        n 

becomes  considerably  more  difficult  to  justify  as  the  number  of  altered 
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terms  increases  and  because  it  was  felt  that  the  additional  complexity 
in  the  resulting  Q  was  self-defeating.   However,  further  research  in 
this  area  may  ultimately  lead  to  a  necessary  and  sufficient  condition 
for  optimality  when  Q  is  not  positive  semidef inite. 

There  remains  but  one  operation  which  may  be  required  to  com- 
plete the  process  of  designing  a  performance  index  for  which  a  given 
closed-loop  pole  configuration  is  optimal.   If  the  state  model  for  the 
plant  (assumed  completely  controllable)  is  not  in  companion  matrix 
canonical  form  it  may  be  desired  to  transform  the  Q  which  has  been 
generated  into  the  coordinate  system  of  the  plant.   This  is  accomplished 
by  the  congruent  transformation, 

T 
Q  =  T  QT, 

where  Q  is  in  the  original  coordinate  system  of  the  plant  and  T  is  the 
non-singular  matrix  of  the  similarity  transformation  (Section  3.4), 

F  =  TFT 

which  transforms  the  original  state  matrix  (F)  into  a  companion 
matrix  (F) . 

It  is  necessary  to  show  that  this  can  be  done  without  loss  of 
generality.   That  is,  that  the  closed-loop  poles  are  invariant  and  that 
no  conjugate  points  are  introduced  under  this  transformation. 
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Theorem  4. 1 

Optimality  is  preserved  under  similarity  transfor- 
mation.  That  is,  if  a  completely  controllable  linear 
plant  and  feedback  control  law  are  optimal  with  respect 
to  a  quadratic  performance  index,  the  system  and  control 
law  subject  to  a  non-singular  coordinate  transformation 
minimizes  the  performance  index  in  the  new  coordinate 
system. 

This  theorem  can  be  easily  proven  by  applying  the  transforma- 
tion required  to  return  the  system  to  its  original  form  from  the 
companion  matrix  formulation  to  the  Riccati  equation.   Then,  by  the 
non-singularity  of  T,  if  the  Riccati  equation  becomes  unbounded  in 
one  coordinate  it  will  in  the  other  as  well. 

4. 3   Design  by  Performance  Index  Iteration 

Once  a  technique,  such  as  the  one  presented  in  the  last  section, 
is  available  for  generation  of  performance  indices  for  optimal  closed- 
loop  pole  configurations  a  simple  design  scheme  becomes  evident.   The 
performance  index  could  be  used  as  a  basis  for  improving  a  first  attempt 
at  design.   In  terms  of  Figure  4.1,  this  would  entail  traversing 
branches  lb  and  2  iteratively,  each  time  altering  the  entries  of  Q  in 
a  manner  intended  to  meet  additional  members  of  the  specification  set. 

This  approach  offers  several  advantages  over  competitive  classical 
design  techniques.   If  the  preliminary  closed-loop  pole  configuration  is 
chosen  so  that 

Y(cjd)  =  | cp  (jcu)}   -  |cp(ja))|   2:  0     for  all  real  cu , 

then  maintenance  of  the  positive  semidef initeness  of  the  resulting  Q 
throughout  all  the  iterations  will  insure  that  the  important  stability 
criteria  of  gain  and  phase  margin  are  at  least  6  db  and  60°,  respectively. 
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This  permits  the  designer  to  concentrate  on  bringing  transient  response 
to  within  desired  limits  and  later  returning  to  the  stability  margins 
if  more  stringent  ones  are  required. 

A  second  advantage  is  that  a  tractable  number  of  design  param- 
eters are  displayed  in  the  Q  matrix.   The  precise  number  can  easily  be 

I  ,  where  the  heavy  brackets  represent  the  greatest 
integer  function;  i.e., 

[a]  =  greatest  integer  £  a. 

This  quantity  does  not  include  the  duplicate  (by  symmetry)  terms  appear- 
ing across  the  diagonal  of  Q. 

This  means  that  for  a  third-order  plant  there  are  at  most  four 
parameters  of  Q  which  are  pertinent  to  the  design;  similarly  for  a 
tenth  order  plant,  at  most  32.   These  numbers  may  seem  rather  large  in 
relation  to  the  number  of  coefficients  of  the  characteristic  poly- 
nomials; however,  these  parameters,  as  will  be  seen  later,  may  be 
related  more  or  less  directly  to  the  transient  response  of  the  closed- 
loop  system,  in  contrast  to  the  classical  design  parameters.   Further, 
far  less  than  the  maximum  number  given  will  be  generally  required;  in 
fact,  usually  it  will  only  be  necessary  to  vary  two  of  them  [H2]. 

It  was  noted  in  the  last  chapter  that  the  first  and  last  terms 

on  the  diagonal  of  equivalent  Q  matrices  must  remain  invariant.   That 

is,  q   must  be  the  same  among  all  members  of  an  equivalence  class  and 

similarly  q   must  not  change.   The  next  theorem  presents  an  important 
nn 

new  result  which  is  a  consequence  of  this  observation. 
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Theorem  4. 2 

If  Q  is  any  member  of  an  equivalence  class  of  weighting 
matrices  (Definition  3.1)  for  a  scalar  nth  order  linear 
system  and  optimal  control  law  in  companion  matrix  canonical 
form,  then 

n    2     n    2 

q„  =  .TT.  R.  -  .TT   r. 
11    1=1   l    1=1   l 

n    2     n    2 
and  q   =  .  E  R.  -  .  £  r.  , 

nn   i=l   i    i=l   i 

where  r    and   R  ,    i=  1,2,3,. ...n 

l        i 

are  the  plant  and  optimal  closed-loop  poles,  respectively. 

Proof 

In  the  discussion  of  the  diagonal  Q  case  in  Section  3.6,  it  was 

2 

shown  that  the  coefficients  of  Y(uu  )  form  the  diagonal  entries.   The 

optimization  problem  implied  by  this  Q  may  have  conjugate  points. 
Nonetheless,  since  the  control  law  of  the  hypothesis  was  specified 
to  be  optimal,  there  must  exist  at  least  one  weighting  matrix  which 
forms  a  performance  index  minimized  by  the  control  law  and  further 
it  must  bea  Y-invariant  of  the  diagonal  Q.   Then,  since  the  first  and 

last  entries  on  the  diagonal  of  invariant  matrices  cannot  change,  these 

2 

entries  must  be  exactly  the  first  and  last  coefficients  of  Y(u)  ). 

Consider  the  monic  polynomial  with  real  coefficients, 

,  ■.     n        n-1  ... 

p(s)  =  s   +  a   „ s    +  ...  +  as  +  a  ,  (4) 

n-1  1     o 

and  its  magnitude  square 

I  /■  xl2     2n   ^    2(n-l)        .   2 
p(juj)    =  uu   +  b   „u)      +  ...  +  b,uu  +  b  .      (5) 
1      '  n-1  1      o 
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By   definition 

Y(u>)    =    |cpk(jw)|2   -    |cp(jw)|2;  (6) 

that  is,  Y  is  the  difference  of  two  polynomials  similar  to  (5). 
Then  the  constant  term  of  Y,  Q11,  is  the  difference  of  two  terms 
corresponding  to  b   in  (5).   Clearly, 

o    o 

and  by  a  well-known  theorem  of  elementary  algebra,  the  constant 

term  of  a    monic    polynomial  is  (to  within  the  sign)  the  product 

of  its  zeros.   Then  b   is  the  square  of  the  product  of  the  roots  of 
o 

(4)  and  application  of  this  abstraction  to  (6)  demonstrates  the 

expression  for  q   . 

The  coefficient  of  cu       of  Y ,  q    ,  is  the  difference  of 

n,n 

terms  similar  to  b     in  (5).   By  examining  the  product  of  p(jua)  and 
its  complex  conjugate  it  is  easily  verified  that 

n-l    n-l     n-2 

Again  it  is  known  from  elementary  algebra  that  for  a  monic  poly- 
nomial such  as  (4) ,  a    is  the  negative  of  the  sum  of  the  zeros  of 
n-l 

4 
(4)  and  a    is  the  sum  of  all  combinations  of  products  of  the  roots 

of  (4)  taken  two  at  a  time.   That  is, 


4 
Here  combination  is  used  in  the  sense  of  combinatorial 

analysis;  i.e. ,  no  distinction  is  made  between  the  ordered  pairs 

(c.d)  and  (d,c). 
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and 


n 
=  -  E  r. 
i=l   ] 


n  i-1 


a^   =   E   E  r.r  .  ,   where  r.  , i= 1,2,3 n  are  the 

n_2   i=2  j=l   X  J 

roots  of  (4)  .   The  expression  for  a   „  can  be  rewritten  by  forming:  the 

n-2  *  b 

sum  of  all  the  possible  permutations  of  products  of  two  roots  (includ- 
ing with  themselves)  and  subtracting  off  the  sum  of  the  roots  squared 
and  dividing  by  two  to  account  for  the  combinations  being  summed  twice; 
i.e., 


1 

*n-2  ~  2 


[(A  o2  -  A  <] 


Then,  using  this  expression  for  a    and  substituting  into  (7)  results 


b   ,  =   E  r2  . 
n-1    .  ,   i 

The  expression  for  q    then  follows  directly  from  application  of  this 
nn 

result  to  equation  (6). 

The  interpretation  of  the  effect  of  q   on  a  system  can  be 

nn       J 

greatly  enhanced  through  the  use  of  a  simple  mechanical  analogy. 

If  the  poles  (r.)  of  a  linear  system  are  thought  of  as  unit  point 

masses,  the  expression, 

n    2 
I  =   E  r2, 

i=l  X 

would  be  the  resulting  moment  of  inertia   with  respect  to  a  perpendic- 
ular axis  of  rotation  through  the  origin   such  a  mechanical  system 

would  possess.   Then  q    represents  the  amount  by  which  the  "moment 
nn 
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of  inertia"  of  the  system  has  increased  with  the  addition  of  optimal 
feedback. 

If  q    is  non-negative  the  optimal  closed-loop  system  will 
nn 

have  at  least  one  pole  further  removed  from  the  origin  than  the  plant 

poles.   That  is,  for   q   >  0   the  closed-loop  system  will  tend  to  be 

*  nn 

faster  in  response  than  the  open-loop,  and,  in  general,  the  larger  q 

is  made  the  faster  the  optimal  system. 

Although  q   does  not  admit  to  analogy  as  well  as  q   ,  its 

effect  on  the  optimal  pole  locations  is  similar  to  that  of  q   j  that 

nn 

is,  the  speed  of  response  also  tends  to  increase  with  q   . 

Recently  there  has  been  a  resurgence  of  experimental  investi- 
gations of  the  effect  of  entries  of  the  Q  matrix  on  the  dynamic  response 
of  the  closed-loop  system  [e.g.,  H2.R1].   Although  it  is  hopeless  to 
attempt  to  glean  a  generalized  design  scheme  from  such  studies,  since 
they  are  bound  to  specific  plant  configurations  (i.e. ,  a  fixed  number 
of  poles  and  zeros) ,  it  is  possible  to  obtain  some  guidelines. 
Houpis  and  Constantinides  [H2]  observed  that  q   and  q   together  have 
primary  effect  on  rise  time  (t  )  and  setting  time  (t  )  and  the  remainder 
of  the  entries  primarily  influence  overshoot.   These  observations  seem 
to  coincide  quite  well  with  the  analysis  given  here. 

At  this  point  a  design  example  would  help  to  solidify  the 
preceding  remarks. 
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Example  4.2 

It  is  desired  to  compensate  the  unstable  plant 


1.0 


s(s-l)(s+2) 


so  that  the  following  performance  specifications  are  met: 

1.  Gain  Margin:   GM  ^  12  db 

2.  Phase  Margin:   cpM  ^  60° 
and  in  response  to  a  step  input, 

3.  Overshoot:   0s  £  5% 

4.  Rise  time  (within  90%):  t   £  1  sec 

r 

5.  Setting  time  (within  1%)  :  t  £  2  sec. 


The  approach  taken  here  will  be  to  develop  a  preliminary 
design,  using  dominant  roots  techniques,  which  is  also  optimal,  and 
to  then  iterate  on  the  entires  of  the  Q  matrix  until  a  suitable  final 
design  is  arrived  at.   By  dealing  only  with  positive  semidefinite  Q' s , 
phase  margin,  and  probably  gain  margin  as  well,  will  be  guaranteed, 
leaving  only  the  dynamic  response  specifications  to  be  investigated. 

Using  the  dominant  root  philosophy,  the  complex  poles  should 
have  a  time  constant  of  approximately  2  seconds  to  meet  the  require- 
ment for  t   and  a  damping  of  approximately  0.7  in  order  to  keep  max- 
imum overshoot  within  5%  and  still  meet  rise  time  specifications 
[D4,  Fig.  4-3,  p.  91].   The  third  pole  removed  by  2.5  times  the  dominant 
time  constant  should  insure  dominance.   The  resulting  preliminary 
design  pole  configuration  is  shown  in  Figure  4.5. 
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2j 


-5   -4   -3    -2   -1 


x  -  plant  poles 

□  -  proposed  closed-loop  poles 


Figure  4.5  Preliminary  Design  Pole  Configuration 


The  algorithm  of  the  preceding  section  is  now  applied  to  obtain 
a  Q  matrix  (if  one  exists).   The  system  is  indeed  optimal  and  a  Q  matrix 
which  will  work  is 


Q  = 


1600  0  0 
0  60  0 
0      0     20 


The  transient  response  of  this  system  was  computed,  using  the  quadratic 
optimization  and  simulation  program  LQL  [B7,B9].   The  important  char- 
acteristics of  the  response  are  recorded  in  the  first  row  of  the  table 
of  Figure  4.6.   It  is  clear  that  response  is  too  slow,  hence,  q   is 
increased;  however,  continued  increase  much  beyond  the  value  of  entry  3 
causes  the  overshoot  specification  to  be  exceeded.   Reducing  q„_  has  the 
effect  of  reducing  the  damping  on  the  dominant  poles  (decreasing  t  )  and 
bringing  the  third  pole  in  toward  the  origin  and  results  in  more  rapid 
settling  (decreasing  t  ).   Run  number  6  meets  the  specifications  on 
transient  response. 
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Run 

t 
r 

Os 

t 
s 

No. 

qil 

q22 

q33 

(sec) 

(%) 

(sec) 

1 

1600 

60 

20 

1.2 

3.3 

2.6 

2 

576 

60 

20 

1.6 

2.0 

3.2 

3  i 

3000 

60 

20 

1.  1 

4.0 

2.3 

4  j 

2500 

60 

10 

1.1 

4.1 

2.2 

5 

2500 

60 

4 

1.  1 

4.0 

2.1 

6  ; 

3000 

60 

4 

1.0 

4.6 

2.0 

7  . 

5000 

60 

4 

0.9 

5.3 

1.8 

Figure  4.6   Transient  Response  vs  q 

li 


The  stability  margins  only  remain  to  be  investigated. 
Since  the  Q  for  all  the  iterations  is  positive  semidef inite, 
GM  ^  6  db  and  cpM  ^  60°.   The  Nyquist  diagram  of  Figure  4.7  shows 
that  the  gain  margin  is  considerably  in  excess  of  12  db;  thus 
iterative  number  6  meets  all  the  design  specifications. 
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Tigurc  -1.7   Xyquist  Plot  of  Final  Design 


Figure  4.8   Root-Locus  Plot  of  Final  Design 

X  -  Plant  poles 

D  -  Closed-loop  poles 

0  -  Feedback  zeros 


90 


Y(t) 


1.0 


Time  (sec) 
Figure  4.9   System  Responses  to  a  Step  Input  -  Preliminary  and  Final  Designs 


91 


The  final  design  meets  the  performance  criteria: 

GM  =  15.18  db 
cpM  =  62.9° 
tr  =   1.0  sec 
t   =   2.0  sec 
•  Os  =   4.6% 

and  has  the  state  feedback  control  law 


54.77 


k  = 


33.53 
7.49 


Although  Figure  4.6  only  tabulates  seven  iterations  on  the  Q 
matrix,  fourteen  iterations  using  program  LQL  were  actually  made, 
requiring  approximately  6  seconds  of  CPU  time  on  an  IBM  360/65  at 
a  cost  of  about  $0. 65. 

This  example  illustrates  two  facets  of  this  method.   First, 
there  are  no  concrete  rules,  as  such,  for  the  manipulation  of  the 
entries  of  Q  to  obtain  desired  changes  in  response,  only  general 
guidelines.   Secondly,  the  speed  and  low  cost  with  which  the  iter- 
ations can  be  computed  make  this  process  of  practical  value. 

A  final  feature  of  this  scheme  is  that  the  resulting  design 
is  not  bound  to  any  specific  form  of  compensation:  this  can  be 
regarded  at  times  either  as  an  advantage  or  a  handicap    However,  with 
the  rapid  advance  being  made  in  automated  compensator  design  [e.g. ,  PI. 
P2] ,  the  disadvantages  associated  with  the  compensator  not  being 
prespecified  are  becoming  largely  illusionary. 
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4. 4   Design  by  Explicit  Performance  Index  Specification 

If  the  design  problem  is  more  fully  specified,  a  useful  varia- 
tion of  the  preceding  scheme  will  often  be  profitable.   For  instance, 
when  specifications  of  the  sort  required  for  Example  4.2  are  given  and 
in  addition  it  is  desired  to  minimize  (in  a  mean-square  sense)  a  spe- 
cific system  output,  or  a  weighted  sum  of  outputs,  the  technique  of 
performance  index  iteration  takes  on  a  particularly  interesting  form. 
Consider  the  scalar  input,  multioutput  linear  system, 

x  =  Fx  +  gu 
T 

y  =  h  x, 

for  which  it  is  desired  to  determine  a  feedback  control  law  so  that 
a  given  (but  at  the  moment  arbitrary)  set  of  classical  performance 
specifications  are  satisfied  and  in  addition  the  performance  index, 


J  =  J    (cyTAy  +  u2)dt,  (8) 


is  minimized,  where  c  is  an  unspecified  scalar  and  A  a  symmetric 
weighting  matrix.   Then  in  the  usual  notation, 

Q  =  cHAHT 

and  cY(id)   =  cY(HAHT;uj)  , 

which  can  be  easily  computed  using  Lemma  3.4.   Any  closed-loop  system 
which  minimizes  (8)  will  have  a  characteristic  polynomial,  cp  (s) , 
such  that 
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l<Pk(Jw)|2  =    |cp(jw)|2  +   cY(u>)  (9) 

by  Theorem  3.3.      Then   if   the  substitutions 


Pk(a)    =    \P^3W)\2 
P(a)   =    |q>(J(u>  | ^ 

and  Y(cr)    =  YGo) 


2    ' 

a=a) 


2    ' 

a=co 


2 

are  made  into  (9) ,  the  resulting  expression, 

Pk(a)  =  p(a)  +  cY(a) , 

is  reminiscent  of  the  classical  root-locus  form,  except  that  the  poles 
of  interest  are  not  those  located  on  the  locus  but  rather  their 
(negative)  square  roots.   Such  plots  have  been  called  "Root-Square- 
Locus"  plots  by  Chang  [C3]  who  made  use  of  them  in  a  different  context 
to  determine  solutions  to  the  so-called  ISE  problems. 

Since  the  root-square-locus  determines  where  the  square  of  the 
closed-loop  poles  may  lie  to  minimize  (8) ,  the  parameter  c  is  the  only 
variable  available  for  use  in  reconciling  any  remaining  performance 
specifications.   If  the  remaining  specifications  include  bandwidth  and 
steady-state  error  (for  a  constant  input)  the  region  of  permissible 
values  for  c  may  be  significantly  reduced. 

If  A  is  positive  semidef inite ,  then  there  exists  a  unity  rank 

Q  such  that 

T    T 

Q  =  hh  -HAH 

T 
and,  further,  h  S(s)  =  \|/(s)  , 


9  1 


is  the  Hurwitz  spectral  factor  of  Y(co) .   Then  the  multioutput  case  can 
be  reduced  to  an  equivalent  single  output  case  if  A  is  positive  semi- 
definite.   That  is,  for  the  purposes  of  the  design,  the  system  can 
be  envisioned  as: 


x  =  (F-gk  )x 
+  gu 


y  =  h  x 


When  y  is  the  equivalent  single  output  derived  above,  the  open-loop 
(k=0)  and  closed-loop  transfer  functions  are,  respectively: 


o  °  cp(jco) 


and 


G  (ju>)  = 


9k(juj) 


It  is  reasonable  to  expect  that  any  bandwidth  or  steady-state 
error  specifications  given  are  in  terms  of  this  equivalent  output;  in 
the  single  output  case  this  would  certainly  be  true.   The  multiple 
output  case  may  require  a  redefinition  of  (8)  to  utilize  the  follow- 
ing procedure,  depending  on  the  interpretation  given  to  bandwidth  and 

steady-state  error  when  more  than  one  output  is  considered. 

5 
Then  the  bandwidth  of  the  closed-loop  system  is  a)   such  that, 

2    1 1 I  2 


K  (Jco  )|   =  i|G  (0>r, 


or  equivalently , 


Y(u)  ) 
o 


2'  c 
1    t(0) 


ivj^l2" 2  K(0)'2 


Clearly  G  (s)  must  be  low-pass,  i.e.  ,  i]j(s)  has  no  zero  at  the 
origin,  for  this  definition  to  be  meaningful. 
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But  |«p    <Ju)>  /      =    |cp<jcu>  |      +   cY(u));    then 

Y(u>   )  .  Y(0) 

o  1 


|cp(ju;o)|2+cY(uJo)      2    |cp(0)|2  +   cY(O) 

Denoting  Y(0)  as  p   and  cp(0)  as  a   (the  coefficients  of  the  least  powers 
in  Y  and  cp)  and  solving  for  c  results  in 

|2         „    2 


|cp(jtOo)|  2ax 


Y(u>   )  p, 

o  1 

P1  ^  0  (see  footnote  5) ,  which  can  be  further  simplified  to 

(10) 


iG(j«,)r    pi 


Then  the  c  in  (8)  necessary  for  the  optimal  system  to  possess  a  speci- 
fied bandwidth,  to  ,  can  be  easily  evaluated  from  expression  (10)  and 
is  a  function  only  of  the  frequency  response  of  the  plant. 

If  bandwidths  much  in  excess  of  the  bandwidth  of  the  plant  are 
required   |G  (jcu  )|  will  become  small,  c  will  rapidly  increase  and  the 
feedback  gains  necessary  to  accommodate  the  resulting  large  coefficients 
of  ^.(s)  will  likewise  increase.   Clearly,  if  the  amount  of  feedback 
gain  is  limited,  the  bandwidth  of  the  optimal  system  cannot  be  much 
larger  than  that  of  the  plant. 

A  similar  procedure  results  in  an  expression  for  steady-state 
error  for  a  unit  step  input.   If  steady-state  error  is  defined  as  the 
difference  between  the  closed-loop  system  input  and  output  as  t  -*  co, 
then  using  the  final  value  theorem  of  transform  theory,  the  final 
value  of  error  (FVE)  is  simply 
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FVE  =  1  -  G  (0) , 


FVE  =  1  -  {— 

ai  +  CP1 


The  square  root  is  allowed  and  has  the  correct  sign  since  i}r(s)  and 
cp  (s)  are  Hurwitz  polynomials.   Then,  by  simple  manipulation,  the 
c  required  for  a  specified  FVE  is 


(1-FVE)2   Pl 


(11) 


for   p   ^  0   and  FVE  ^  1. 

As  should  be  expected,  steady-state  error  and  bandwidth  are 
conflicting  requirements.   Whereas  large  bandwidths  result  in  large 
values  for  c,  increased  accuracy  follows  from  decreasing  c.   If  zero 
FVE  is  specified,  then  from  (11) 

c  =  1  -  &1/p   . 

Substituting  this  value  for  c  into  (10) 

2 

1 


|o  (jo.  )|2      Pl 

1   O     O  ' 


If  the  required  bandwidth,  uu  ,  is  greater  than  the  open-loop  bandwidth, 

o 

then 

P 


lGo(jV'2<i4 
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and  |G0(jVI      =   2<72 

Pl+31    2ai 

hence,  it  is  possible  to  obtain  closed-loop  bandwidth  greater  than  that 

of  the  plant  and  zero  FVE  simultaneously  if 

.   2 
Pl  >  31  • 

This  condition  can  always  be  met  by  providing  additional  gain  in  the 
plant  path.   That  is,  if  a  is  the  constant  coefficient  of  the  numer- 
ator of  G  (iu))  and  g  is  a  variable  gain  on  the  output  of  the  plant, 
o 

then  g  is  simply  increased  until 

2    2 
Pl  "  (gao)   >  31  * 

Although  the  preceding  treatment  has  dealt  exclusively  with 

a  unity  rank  Q,  the  actual  computation  of  the  required  optimal  control 

T 
law  may,  of  course,  be  carried  out   using  any  Q  equivalent  to  chh  . 

These  methods  for  computing  optimal  systems  satisfying  steady-state 
error  and  bandwidth  specifications  may  also  be  profitably  used  as  an 
initial  design  for  the  performance  index  iteration  techniques 
(Section  4.3)  when  additional  criteria  are  to  be  met. 

In  summary,  when  it  is  desired  to  remove  the  disturbances  from 
a  specific  system  output  or  weighted  sum  of  outputs,  additional  tech- 
niques become  available  for  selection  of  the  correct  performance  index 
to  meet  design  requirements.   The  squares  of  the  permitted  closed-loop 
locations  may  be  identified  by  the  use  of  the  root-square-locus,  and 
a  performance  index  generated  to  place  them  anywhere  along  that  locus. 
If  steady-state  error  and/or  bandwidth  specifications  are  to  be  met, 
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performance  indices  may  be  simply  constructed  which  result  in  optimal 
systems  satisfying  these  criteria.   In  each  case,  these  performance 
indices  are  explicitly  defined. 

4. 5   Sampled-Data  Controller  Design 

This  section  will  consider  a  design  problem  which  is  particularly 
well  adapted  to  solution  by  techniques  arising  from  the  inverse  problem. 
The  design  of  digital  controllers  for  continuous  systems  is  attacked, 
using  methods  which  represent  a  new  approach  to  a  long-standing  problem. 

The  specific  problem  to  be  dealt  with  is  identical  to  the  one 
considered  earlier  in  this  chapter  with  the  exception  that  the  result- 
ing compensator  is  constrained  to  be  a  digital  device  rather  than 
continuous.   That  is,  it  is  desired  to  design  a  digital  compensator  for 
the  linear,  completely  controllable,  scalar  input  system 

x  =  Fx  +  gu  (12) 

such  that  a  given  set  of  classical  performance  specifications  are 
satisfied. 

The  approach  will  be  to  use  the  techniques  of  the  preceding 
sections  to  design  a  performance  index,  having  a  corresponding 
continuous  optimal  control  law  which  meets  the  performance  specifi- 
cations.  The  sampled-data  model  for  the  continuous  system, 

x   .  =  Ax  +  Iu  ,  (13) 

m+1     m     m 

is  then  computed  so  that  the  solutions  to  (12)  and  (13)  agree  at  the 
sample  points  when  the  input  is  held  constant  over  the  sample  interval, 
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as  it  would  be  if  it  were  the  output  of  a  zero-order  hold  of  a  digital 
controller.   The  required  matrices  A  and  F  in  (13)  can  be  determined 
[B7]  from  the  solution  to  (12),  i.e., 

t 

P 

t 


X(t)  =  §(t-t  )x(t  )  +  P   §(t-T)gu(T)  dT, 
o     o     \ 


Ft 
where  $(t)  =  e   , 

by  evaluating  it  over  the  interval,  (t  =  mAT,  t  =  (m+l)AT),  where 
AT  is  the  sample  period.   Then  with  the  convention  that  variables  eval- 
uated at  time  mAT  be  denoted  by  a  subscript  m, 

A  =  $(AT) ,  (14) 

(m+1)  AT 
and  r  =  J        §((m+l)AT-T)dTg.  (15) 

mAT 

If  §  in  (15)  is  written  as  the  matrix  exponential  and  expanded  in  a 
Maclaurin's  series,  the  integral  can  readily  be  evaluated  as  the 
infinite  series, 

r - I  m^  **  • 

The  expressions  for  A  and  T,    (14)  and  (16),  are  ideal  for  direct 
numerical  evaluation  when  (14)  is  expanded  in  a  Maclaurin's  series. 

The  Q  derived  from  the  continuous  design  earlier  is  then 
inserted  in  the  discrete  performance  index, 


J  =   £   (xTQx.+u2) ,  (17) 

j=0   J  J  J 
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and  the  optimal  control  law  for  the  sampled-data  model  (13)  and  the 
discrete  performance  index  (17)  is  computed.   This  optimization  problem 
|"K3]  is  the  discrete  analog  of  the  continuous  one  considered  throughout 
this  work.   This  problem  results  in  a  constant  optimal  feedback  law 
computed  similarly  to  the  continuous  case.   If  the  optimal  sampled-data 
control  law  is, 

u   .,  =  -  k  x  , 

m+1        m 

then  it  is  hoped  that  the  closed-loop  sampled-data  system  will  perform 
similarly  to  the  continuous  closed-loop  system  which  was  originally 
designed,  so  that  the  specifications  are  complied  with  by  the  discrete 
design  as  well. 

As   AT  -*  0   it  is  clear  that  the  optimal  discrete  system  will 
respond  identically  to  the  optimal  continuous  design.   For  AT  suffi- 
ciently small  the  discrete  design  should  perform  within  the  specifica- 
tions.  If  the  sample  rate  is  too  slow  to  faithfully  reproduce  the 
desired  response,  the  Q  generated  from  the  continuous  design  will 
nonetheless  provide  a  basis  for  design  by  performance  index  iteration, 
as  was  done  for  the  continuous  case  in  Section  4.3. 

By  way  of  illustration,  the  example  of  design  by  performance 
index  iteration  will  be  repeated,  this  time  with  the  goal  of  designing 
a  digital  compensator  which  meets  the  prescribed  time  domain  specifi- 
cations.  The  frequency  domain  specifications  (gain  and  phase  margin) 
will  not  be  investigated  since  these  terms  require  redefinition  for 
sampled-data  systems.   Frequency  domain  analysis  for  sampled-data 
systems  can  generally  be  accommodated  with  no  difficulty  through  the 
use  of  the  "bilinear  transforms"  [K12]. 
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Example  4. 3 

The  plant  of  Example  4.2  is  to  be  compensated  through  the  use 
of  a  sampled-data  controller,  with  an  (at  present)  unspecified  sample 
rate,  so  that  it  meets  the  following  time  domain  performance  cri  ia. 
In  response  to  a  step  input  the  system  should  have 

1.  Overshoot:   0s  £  5% 

2.  Rise  time  (within  90%):   t  £  1  sec 

r 

3.  Settling  time  (within  1%) :   t  £  2  sec. 
In  the  earlier  example 


Q  = 


3000 

0 

0 

0 

60 

0 

0 

0 

4 

was  found  to  result  in  an  optimal  control  law  which  met  these  specifi- 
cations for  the  system  in  companion  matrix  canonical  form. 

All  that  remains  is  to  compute  the  discrete  model  of  the  system 
and  optimal  discrete  control  law  corresponding  to  the  Q  above.   This 
was  done,  again  using  program  LQL,  with  sample  period  of  0.05,  0.1, 
and  0.2  seconds.   All  three  sample  rates  resulted  in  compensated 
systems  which  easily  met  the  specifications.   The  step  response  of 
the  optimal  discrete  system  with   AT  =  0.1   is  representative  of  the 
three  trials.   This  response  was  within  approximately  1%  (of  the  step 
magnitude)  of  the  response  of  the  final  design  in  Example  4.2  for 
the  first  1.5  seconds  and  within  less  than  0.1%  thereafter. 
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The  three  sample  rates  were  chosen  to  be  near  the  "rule-of- 
thumb"  value  of  approximately  one-tenth  the  fastest  time  constant  of 
the  optimal  continuous  system  and  are  therefore  reasonable  values. 
A  sampling  period  of  one-half  second  was  also  used  and  resulted  in 
a  closed-loop  system  response  very  similar  to  that  desired  but  just 
failing  to  meet  the  specifications. 

It  would,  of  course,  be  rather  extravagant  to  employ  a  digital 
compensator  to  simply  feed  back  a  linear  combination  of  the  states. 
However,  with  the  availability  of  an  arithmetic  unit  (and  presumably 
some  memory)  in  the  compensator,  the  implementation  of  a  Kalman  filter 
to  estimate  the  states  in  the  presence  of  noise  is  a  relatively  minor 
extension. 

As  mentioned  earlier,  the  original  motivation  for  this  entire 
study  was  to  develop  a  technique  for  retrofitting  a  digital  compensator 
to  a  system  with  a  satisfactory  continuous  compensator  without  altering 
system  performance.   It  is  now  clear  how  this  can  be  accomplished. 
The  techniques  of  Section  4.2  are  applied  to  the  original  continuous 

system  to  determine  the  performance  index,  if  any,  minimized  by  the 

7 
original  configuration.    If  the  original  system  is  optimal ,  then  the 

design  of  a  digital  compensator  proceeds  as  in  the  example.  When  sign 

indefinite  Q' s  are  admitted,  this  is  not  a  severe  restriction. 


7 
Compensator  poles  should  be  removed  from  both  open-  and  closed- 
loop  systems  before  optimality  is  tested  or  the  performance  index  con- 
structed. 


CHAPTER  V 


CONCLUSIONS 


5. 1   Summary  of  Results 

It  often  seems  in  the  pursuit  of  a  solution  to  a  specific 
problem  that  the  many  incidental  and  preliminary  results  one  encounters 
along  the  way  are  of  as  much  lasting  value  as  those  originally  sought. 
That  apparently  has  been  the  case  in  this  study. 

The  original  intent  of  the  research  presented  here  was  to 
develop  a  scheme  to  use  optimal  control  theory  for  the  automated  design 
of  digital  compensators  to  replace  existing  continuous  ones.   This 
leads  to  an  investigation  of  the  inverse  optimal  control  problem,  to 
determine  when  a  solution  to  the  inverse  problem  exists  and  how  it 
is  obtained.   It  was  found  that  the  criterion  for  optimality  is  a  fre- 
quency domain  result  and  leads  naturally  to  consideration  of  the  mean- 
ing of  optimality  in  a  classical  (time  and  frequency  domain)  context. 
From  the  study  of  the  common  ground  of  classical  and  optimal  control 
theory,  the  possibility  arose  of  using  the  linear  quadratic  loss 
problem  as  a  computational  vehicle  for  design  of  linear  controllers 
(continuous  and  discrete)  to  classical  performance  specifications. 

The  contributions  of  this  research  can  be  summarized  as  follows: 

1.  A  complete  theoretical  investigation  of  the  inverse 
optimal  linear  regulator  problem  for  quadratic  per- 
formance  indices  with  positive  semidefinite 
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state  weighting  matrices  and  scalar  input  systems  is 
reviewed  and  preliminary  results  for  the  case  where  the 
weighting  matrix  is  allowed  to  be  sign  indefinite  are 
given. 

2.  The  equivalence  of  performance  indices  for  scalar  input 
linear  quadratic  loss  problems  is  resolved  and  a  proce- 
dure for  generating  the  entire  equivalence  class  of  cost 
functions  equivalent  to  a  given  performance  index  is 
developed. 

3.  Practical  numerical  methods  for  determination  of  system 
optimality  and  computation  of  solutions  to  the  scalar 
input  inverse  problem  are  discussed. 

4.  Some  of  the  effects  of  specific  elements  of  the  perform- 
ance index  on  optimal  system  performance  and  pole  loca- 
tions are  determined. 

5.  The  problem  of  designing  optimal  systems  to  meet  classical 
performance  specifications  is  encountered  and  some  defin- 
itive results  obtained. 

6.  A  solution  to  the  problem  of  designing  a  sampled-data 
controller  to  approximate  the  performance  of  continuous 
controller  is  given. 

5.2   Suggestions  for  Future  Research 

Many  of  the  results  of  this  work  are  quite  open-ended.   For 
instance,  the  goal  of  using  optimal  control  theory  for  completely 
automated  design  of  control  systems  subject  to  classical  performance 
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specifications  is  still  far  off.   Some  areas  where  research  can  be 
invested  with  immediate  benefit  are  outlined  below. 

1.  Although  some  investigations  of  the  inverse  problem  fbr' 
multi-input  systems  have  been  made,  they  are  generally 
very  cursory  and  much  work  in  the  spirit  of  Chapters  III 
and  IV  here  remains  to  be  done. 

2.  Necessary  and  sufficient  algebraic  conditions  for  the 
non-existence  of  conjugate  points  when  the  state  weight- 
ing matrix  of  the  quadratic  performance  index  is  sign 
indefinite  are  eagerly  awaited. 

3.  More  extensive  relations  between  performance  index  speci- 
fication  and  optimal  system  response  need  to  be  developed. 

4.  Techniques  similar  to  those  of  Chapter  IV  can  probably 
be  profitably  applied  to  the  so-called  "model  following 
problem"  [Al]. 


APPENDICES 


APPENDIX  A 


ADDITIONAL  PERFORMANCE  INDICES  ACCOMMODATED 
BY  THE  THEORY 


A. 1   Introduction 

As  mentioned  in  the  text,  several  performance  indices  other 
than  the  one  specifically  discussed,  i.e., 

J  =  J      (xTQx  +  uTRu)dt,  (1) 

t 
o 

may  be  formulated  which  lend  themselves  well  to  the  interpretations 
and  techniques  of  this  work.   To  demonstrate  this,  four  quadratic 
performance  indices,  which  may  be  viewed  as  variations  of  (1),  are 
shown  to  be  reducible  to  (1)  by  some  rather  simple  transformations. 

None  of  these  performance  measures  actually  extend  any  of  the 
results  presented,  except  perhaps  heuristically.   Their  value  is  that 
they  allow  the  designer  to  initially  specify  a  performance  index  which 
may  be  more  appealing  than  (1)  for  one  reason  or  another  in  a  particular 
case  and  then  to  proceed  to  reduce  it  to  the  form  of  (1)  for  applica- 
tion of  the  design  techniques  presented  here. 

A.  2   Quadratic  Performance  Index 

with  Exponential  Weighting 

An  interesting  variation  of  (1)  discussed  by  Anderson  and  Moore 
[A4]  is  to  include  an  exponential  weighting  term  in  the  integrand  of 
the  functional  (1).   That  is, 
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I*        2ctt      T  T 

J   =  J       e         (xQx+u  Ru)dt,      R   positive  definite  (2) 

t 


and  the  completely  controllable  system  to  be  optimized  is 

x  =  Fx  +  Gu  .  (3) 

If  the  resulting  control  law  and  the  system  are  time  invariant,  then 

it  is  reasonable  to  assume  that  the  states  must  decay  at  least  as 

-cyt 
fast  as  e    in  order  for  (2)  to  remain  bounded.   This  will  indeed  be 

shown  to  be  the  case  for  infinite  final  time;  that  is,  the  real  parts 

of  all  the  closed-loop  poles  must  be  less  than  -a. 

Let 

a    ort      a    at 
x  =  e  x  and  u  =  e  u, 

*    Cft.     0/t 
then  x  =  e  x  +  ae     x, 

reapplication  of  the  definition  of  x  and  substitution  from  (3)  for  x 

result   in 

*•  0/t   ,  A 

x  =    e      (Fx  +   Gu)    +   c/x 
or  x  =    (F  +   o/I)x  +    Gu    .  (4) 

Then  (2)  can  be  rewritten  as 

J  =  /   (xTQx  +  u  Ru)dt 

t 
o 

which  along  with  (4)  is  an  optimal  regulator  problem  in  the  usual  sense 

and  has  the  optimal  control  law  (assuming  no  conjugate  points) , 

u  =  -  KT(t)x. 
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Returning  the  control  law  to  the  original  coordinate  system, 

-at  -T_,  at 
u  =  -  e    K  (t)e  x 

=  -  K  (t)x  , 

and  if  t   is  allowed  to  become  infinite,  K  becomes  a  constant  control. 

If  the  optimal  system  is  stable,  x  approaches  zero  as  t  becomes  infin- 

-Ot 
ite;  hence,  x  decays  at  least  as  fast  as  e    and  the  original  inter- 
pretation of  (2)  is  shown  to  be  correct. 

In  summary,  the  performance  index  with  exponential  weighting 
can  be  accommodated  by  a  redefinition  of  the  plant  which  has  the  effect 
of  shifting  the  plant  poles  by  -a.      If  Ot   is  positive,  this  can  be 
thought  of  as  giving  the  optimization  problem,  which  tends  to  generate 
closed-loop  poles  to  the  left  of  the  plant  poles,  a  "head  start"  of  ex. 

A. 3   Quadratic  Performance  Index 
with  Cross-Product- 

Consider  now  the  generalization  of  (1)  which  includes  cross- 
products  of  the  states  and  inputs,  i.e.  , 

J  =  J        (xTQx  +  2xTAu  +  uTRu)dt,  (5) 

t 
o 

where  A  is  an  n  xm  real  matrix,  m  is  the  number  of  inputs,  and  R  is 
positive  definite  and  symmetric. 

As  before,  a  redefinition  of  the  model  will  reduce  this  varia- 
tion to  the  form  of  (1)  ,  in  this  case  only  the  control  need  be  modified. 
Let 

u  =  u  +  B  x, 


110 


where  B  is  an  n  x m  real  matrix,  then  the  integrand  of  J  is 

T        T   a   T  T   T      T 

x  Qx  +  2x  A(u-B  x)  +  (u-B  x)  R(u-B  x) 

=  x  (Q-2AB  +BRBT)x  +  2xT(A+BR)u  +  uTRu. 
The  problem  will  reduce  to  the  usual  form  if 

A  =  BR; 
since  R  is  positive  definite  it  is  also  non-singular  and  the  required 
B  is 

B  =  AR_1  .  (6) 

The  optimization  problem  with  performance  index  (5)  may- 
be solved  by  replacing  the  Q  of  the  usual  performance  index  (1)  with 

Q  =  Q  -  2ABT  +  BRBT 

from  (5);  substitution  of  (6)  for  B  simplifies  this  considerably  to 

-1  T 

Q  =  Q  -  AR  A   . 

Then,  if  the  customary  optimization  problem  with  Q  has  no  conjugate 

~T 
points,  it  will  have  the  optimal  control  law,  u  =  -  K  x,  and  the  optimal 

A       —IT 

control  for  (5)  is       . u  =  -  (K  +  AR   )  x.  (7) 

Clearly,  there  is  no  increase  in  generality  if  A  of  (5)  is 
chosen  to  be  non-null,  since  no  additional  optimal  systems  would  be 
admitted.   Formulation  (5)  permits  a  change  in  the  appearance  of  the 
optimization  problem  with  no  effect  on  the  theory. 
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A. 4   Quadratic  Performance  Indices 
with  Derivative  Weighting 

Two  variations  of  derivative  weighting  performance  indices  will 

now  be  considered.   First,  a  quadratic  form  composed  of  the  (first) 

derivatives  of  the  states  will  be  added  to  (1);  that  is, 

J  =  J      (xTQx  +  xTCx  +  u  Ru)dt,  (8) 

t 
o 

where  C  is  a  symmetric  n  xn  real  matrix.   Upon  substitution  from  (3) 

for  x  the  integrand  becomes 

T  T  T 

x  Qx  +  (Fx+Gu)  C  (Fx+Gu)  +  u  Ru 

and  simplifies  to 

x  (Q+FTCF)x  +  2xTFTC  Gu+uT(R+GTCG)u. 

This  integrand  is  of  the  same  form  as  the  performance  index  of  the 
last  section  and  hence  can  possibly  be  put  in  the  customary  form  by 
a  simple  redefinition  of  u.   In  this  case  the  required  B  (6)  is 

B  =  FTCG(R+GTCG)-1  (9) 

and  may  not  exist.   If  the  necessary  transformation  exists 

R  +  GTCG 

must  be  non-singular.   This  must  also  be  non-singular  for  the  equiva- 
lent R  of  the  resulting  performance  index  to  be  positive  definite; 
thus,  if  the  optimization  problem  with  (8)  has  a  solution,  the  trans- 
formation (9)  must  exist.   Then  (8)  is  again  different  from  (1)  in 
appearance  only. 
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The  second  variation  of  derivative  weighting  is  to  include 
a  quadratic  form  of  the  derivative  of  the  control;  i.e., 

J  =  J   (x  QX  +  u  Du  +  u  Ru)dt, 
t 


(10) 


where  D  is  a  real  m  xm  symmetric  matrix.   This  case  is  handled  by 
renaming  the  control  and  adding  m  (one  per  input)  additional  states. 
Let  u  =  ii,  then  graphically  the  original  configuration 


Plant 

1 

2 

x 

u 

m 

is  changed  so  that  the  system  to  be  optimized  is 


u1 LVsJ 


n+1 


1/sj 


1/s 


n+2 


Plant 


Figure  A. 1   Modified  Plant 
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Define  the  state  model  for  the  modified  (n  +  m  order)  plant  as 


x  =  Fx  +  Gu, 


where 


F  = 


F   G 
0   0 


(11) 


and  x  has  components  as  in  Figure  A.l.   Then,  if  Q  and  R  in 

J  =  J      (xTQx  +  uTRu)dt 

t 
o 

are  chosen  as 


(12) 


Q  = 


Q   0 
0    R 


and   R  =  D, 


the  solution  to  the  optimization  problem  with  derivative  weighting  (10) 
may  be  obtained  from  the  customary  problem.   Assume  that  the  optimiza- 
tion problem  with  (11)  and  (12)  has  a  solution  (no  conjugate  points) 
which  is 

u  =  -  K  x  .  (13) 

By  partitioning  K  into  two  matrices, 


Kl 

— 

K2 

V 

and  using  the  definition  for  x  from  (11) ,  the  control  law  (13)  can  be 

rewritten  as 

t 

U  =  -  J  '   (K  x  +  K  u)dt. 

t 
o 
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In  this  form  a  realization  for  the  optimal  control  law  in  terms  of  the 
original  plant  is  clearly: 


Controller 


Plant 


■9-  n  variables 
~— -  m  variables 

Figure  A. 2   Synthesis  of  Optimal  Control  for  Performance 
Index  with  Control  Derivative  Weighting 


This  case  is  somewhat  different  from  the  other  performance 
indices  which  have  been  considered  in  that  the  optimal  control  must 
be  synthesized  with  a  dynamic  compensator.   However,  it  is  obvious 
from  the  modified  form  of  the  plant  and  the  performance  index  that 
it  also  is  but  a  variation  of  the  original  problem. 


APPENDIX  B 


PROOF  OF  SUFFICIENCY  OF  THEOREM  3.1 


To  prove  the  sufficient  part  of  Theorem  3.1  it  must  be  shown 
that  for  a  scalar  input  completely  controllable  system, 

x  =  Fx  +  gu, 
and  stable  completely  observable  control  law, 


u  =  -  k  x; 


if 


ll+kT$(jto)g|2  ^   i 


Hs)    =  (sI-F)   , 


for  all  real  uu ,  then  the  control  law  (2)  is  optimal. 

Since,  by  the  hypothesis,  the  system  (1)  is  completely 
controllable  it  can  be  taken  without  loss  of  generality  to  be  in 
companion  matrix  canonical  form  (Section  3.4).   Then  by  Lemma  3.1, 


r~.        *~\ 


S(s)  =  cp(s)§(s)g  = 


n-1 


,  where  cp(s)  =  det  (sI-F) 


(1) 


(2) 


(3) 


(4) 


and 


T        kl  +  k2S  + 
kT§(s)g  =  — ^ 


+  k  s 


n-1 


cp(s) 


K(s) 

cp(s) 


Application  of  this  expression  to  (3)  results  in, 

(2 


|l+k  §(j(ju)g| 


cp(jqj)  +  K(,jxO  , 
cpCjw) 


(5) 
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then  from  the  assumption  of  companion  matrix  canonical  form  for  (1) , 
it  follows  that 

cp(jcu)  +  K(ju)) 

is  the  closed-loop  characteristic  polynomial,  which  will  be  denoted 


9k(ju>) 


Inequality    (5)    can  now  be  rewritten    as 

2 


Y(u))    =    |<Pk<J«>>|       "    |t?(ju))|      *  0, 


displaying  the  non-negativity  of  Y  which  allows  it  to  be  factored 
(Lemma  3.2)  as 

Y(u))  =  i{K-ju))\|f(j(u)  , 

where  i|f(s)  is  a  real  Hurwitz  polynomial  and  can  be  represented  as  the 

vector  product 

T        T 
\!;(s)  =  p  S(s),   p   =  [p  ,p    ..  ,p  ]  . 
1   <£      n 

From  the  proof  of  the  necessary  part  of  Theorem  3.1,  if  k  is  an 
optimal  control  law,  then 

IT       1 2         T  T 
1+k  f(j(ju)g|   =  1  +  g  §  (-juj)Q$(juu)g, 

which  can  be  rewritten  in  the  same  manner  as  (5)  to  result  in 

Y(cju)   =    |cpk(jcu)|2   -    JcpCjoj)  | 2  =  gT§T(-juj)Q$(juj)g.  (6) 

Again  evoking  relation  (4) ,  it  is  clear  that 

Q  =  PPT  (7) 

will  satisfy  (6). 

If  it  can  now  be  shown  that  the  use  of  this  weighting  matrix 
in  an  optimization  with  (1)  will  result  in  k,  the  theorem  is  proved. 
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Suppose  that  the  performance  index, 

00 

f»   T   T      2 
J  =  J  (x  pp  x  +  u  )dt , 

o 

subject  to  the  system  (1)  is  minimized  by  the  control  law,  k  ^k,  and 

denote  the  resulting  closed-loop  characteristic  polynomial  as  cp^1 ( s) . 
By  (6)  and  the  definition  of  Q  in  (7) 

2    i  ,.    .  |2 


K/jw)!  =  |cpkUw)|  ;  (8) 


then  if  k  is  shown  to  be  a  stable  control  law,  the  unique  Hurwitz 
spectral  factor  (Lemma  3.2)  of  both  sides  of  (8)  will  be  identical 
and  k  =  k. 

The  optimal  control  law  k  will  be  stable  if  the  pair  [F,p] 
is  completely  observable  (Theorem  2.5).   Suppose  they  are  not  com- 
pletely observable.   Then  the  numerator  and  denominator  of 

Ts/  ^     Us) 
P  Hs)s   =9Ts7 

will  have  a  common  factor  [K5]  which  must  have  zeros  with  negative 
real  parts  since  ty  is  Hurwitz.   By  rewriting  (6)  as 

|cpk(Jw)|   -  U(j«))|   +  |cpCJui)  |  , 

it  is  clear  that  cp,  (s)  must  contain  the  same  common  factor,  since  k 
k 

is  a  stable  control  law  and  cp  (s)  must  have  only  zeros  with  negative 

real  parts.   Then  from  the  observation  in  (5)  that 

cpk(s)  =  cp(s)  +  K(s), 

K  must  also  contain  the  same  common  factor  as  cp  and  cp  .   But  this 

k 

contradicts  the  stipulation  in  the  hypothesis  that  k  be  completely 
observable.   Hence  the  pair  [F,p]  are  completely  observable  and  the 
theorem  is  proven. 


APPENDIX  C 
THE  INVERSE  PROBLEM  -  NUMERICAL  DETAILS 

C.  1   Introduction 

In  this  appendix  the  numerical  methods  developed  during  the 
course  of  this  research  for  computation  of  solutions  to  the  inverse 
problem  are  discussed.  Specifically,  the  algorithms  presented  here 
deal  with  the  problems  of: 

1.  computing  Y(io)  for  a  given  system  configuration, 

2.  determining  the  non-negativity  of  Y,  and 

3.  the  computation  of  the  Hurwitz  spectral  factor  of  Y 
once  it  has  been  determined  to  be  non-negative. 

These  are  essentially  the  procedures  required  to  determine 
if  a  completely  controllable  system  and  stable  scalar  control  law  are 
optimal  with  respect  to  a  quadratic  performance  index  with  Q  positive 
semidefinite  and  when  it  has  been  ascertained  to  be  optimal,  computing 
a  unity  rank  Q  which  forms  a  performance  index  minimized  by  the  optimal 
system.   The  additional  numerical  operations  required  to  investigate 
the  optimality  of  a  closed-loop  system  with  respect  to  a  performance 
index  which  has  a  Q  not  positive  semidefinite  are  very  straightforward 
(Section  4.2)  and  nothing  would  be  gained  by  including  them  here. 

C.2   Computation  of  Magnitude  -  Square  Polynomials 

The  first  step  in  the  determination  of  the  optimality  of 
a  given  system  configuration  is  computation  of  the  real  even  polynomial , 
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Y(uj)    =    |cp    <jou>  |      -   JcpCjcu)  |      . 

This  necessitates  the  calculation  of  the  coefficients  of  the  magnitude- 
square  polynomials  for  the  open-  and  closed-loop  characteristic  poly- 
nomials; the  algorithm  presented  in  this  section  is  an  efficient  tech- 
nique for  the  computation  of  these  coefficients.   The  method  presented 
is  unique  in  that  it  does  not  require  complex  arithmetic  and  computes 
only  the  non-zero  terms. 

Identification  MAGSQ 

SUBROUTINE  MAGSQ  (X,  Z,  NDIM) 
Purpose 

To  compute  the  magnitude-square  polynomial 

Z(U)2)  =  |x(ju))|2  =  X(ju>)  X(-juu) 

of  a  given  real  polynomial. 

Method 

2 
Let   Z(iu  )  =  X(juj)  X(-jiu) 

N  N 

and  Z(ou2)  =   E  z.u)  X        X(s)  =   E  x  s1 , 

i=0   X  i=0 

2 
then  since  Z(a>  )  is  even  only  the  even  power  coefficients  need  be 

N   N 
computed,  that  is,   z.  ,  =   S   £   (j)   x.(-j)   x 
1+k   i=0  fc=0       x      ^ 

such  that   i  +  k   is  even.   Each  coefficient  of  Z  is  then, 


k  .  . . i+k 


z.   =       (-1)"  (j)*"~  xiXk 


i+k=£ 

even 

and  because  i  +  j  is  an  even  integer, 

...i+k    .  lN(i+k)/2 
(j)    =  (-1) 
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and  the  previous  expression  may  be  rewritten 

z  .  =        E   (-1)   (-1)        x.x   . 
*   .  ,  .  l  k 

i+k=£ 

even 

Then  the  expression  used  for  computing  the  coefficients  of  Z  is 

„    ,  1,(i+3k)/2 
z  .  =        S   (-1)         x.x  , 

i+k=X 

even 

which  requires  no  complex  arithmetic  and  only  the  N+l  non-zero  terms 

are  computed. 

Usage 

CALL  MAGSQ  (X,  Z,  NDIM) 
Input 

X  -  vector  of  polynomial  coefficients  stored  from  smallest  to 

largest  power  terms  (single  or  double  precision,  see  Remarks) 
NDIM  -  number  of  coefficients  of  X  and  Z,  i.e.  ,  N+l 
Output 

Z  -  vector  of  magr  tude-square  polynomial  coefficients  stored 
from  smallest  to  largest  even  power  terms,  i.e.  ,  nth  term 
of  vector  is  coefficient  of  uu       (single  or  double  pre- 
cision, see  Remarks) 
Remarks 

Coefficient  arrays  X  and  Z  are  presently  double  precision  but 
may  be  made  single  precision  by  placing  a  "c"  in  column  one  of  card  12 
as  indicated  below: 

C    DOUBLE  PRECISION  X,Z 
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C.3   Test  for  Non-negntivity  of  Real, 
Even  Polynomials 

After  the  computation  of  the  magnitude-square  of  the  open- 

and  closed-loop  characteristic  polynomials,  optimality  of  the  system 

with  respect  to  a  quadratic  performance  index  with  positive  semidefi- 

nite  state  weighting  matrix  can  be  determined  from  the  sign  definite- 

ness  of  Y(uj)  .   This  test  is  provided  by  subprogram  NONNEG  and  supporting 

subroutine  NSIGNV;  both  subroutines  are  Fortran  integer  functions. 

The  heart  of  the  procedure  is  a  novel  method  for  computation  of  a 

normalized  Routh  Array  and  the  only  known  implementation  of  Siljak's 

test  for  non-negativity. 

Identification  NONNEG 

FUNCTION  NONNEG  (P ,NDIM, IPRT) 
Purpose 

To  verify  non-negativity  of  even  real  polynomials,  that  is, 
to  determine  the  optimality  with  respect  to  a  quadratic  performance 
index  (with  positive  semidefinite  state  weighting  matrix)  of  a  constant 
coefficient  linear  dynamical  system  with  a  given  scalar  feedback 
control  law. 
Method 

The  requirement  that  the  real  even  polynomial,  Y(u;)  ,  be  non- 
negative  for  all  real  to  is  equivalent  (Section  4.2)  to  requiring  that 

2 

Y(u)  )  possess  no  positive  real  zeros  of  odd  multiplicity.   Hence  this 

algorithm  relies  on  properties  related  to  root  location  to  determine 
non-negativity. 

The  procedure  begins  with  four  simple  necessary  tests  (one 

sufficiency  test  is  made  incidental  to  a  necessary  test) ;  if  all  are 
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successful  (and  the  sufficiency  test  fails)  a  necessary  and  sufficient 
condition  for  positive  definiteness  is  tested.   If  this  test  fails, 
a  final  necessary  and  sufficient  test  for  non-negativity  is  made. 
In  general,  the  philosophy  of  this  procedure  is  to  attempt  to  resolve 
the  sign  definiteness  of  the  given  polynomial  with  the  simplest  test 
possible  and  to  resort  to  sophisticated  techniques  only  if  all  else 
fails.   The  following  necessary  conditions  are  first  tested  (in  order): 

i=l  x 

1.  a  £  0,     i.e. ,     Y(0)  £  0 

2.  axT  ,  £  0,   i.e.  ,      lim  Y(uj)  a  0 

N+l 

(l)-»00 

2 

3.  Y(cu  )  must  have  an  even  number  of  sign  variations  between  succes- 
sive terms  (a.'s)  in  order  to  have  an  even  number  of  positive 

l 

real  zeros  by  Decartes1  rule  of  signs.   (If  there  are  zero  sign 
variations  the  coefficients  are  all  non-negative  and  this  con- 
dition is  also  sufficient. )   The  test  for  an  even  number  of  sign 
variations  is  actually  implicit  in  1. and  2  (above);  however,  it 
is  made  as  a  convenient  test  for  zero  sign  variations  and  as 
a  check  for  possible  numerical  problems. 

If  all  the  preceding  tests  on  necessary  conditions  are  passed, 
a  normalized  Routh  Table  TG4]  is  constructed  by  normalizing  the  rows 
to  have  a  maximum  modulus  of  unity  and  the  following  sufficiency  test 

[J2]  is  made: 

2 
A  real,  even  polynomial  Y(w  )  has  no  positive  real  zeros  (hence 

is  positive  for  all  real  cu)  if  and  only  if  the  number  of  sign  variations 
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between  successive  coefficients  is  even  (3  above)  and  the  number  of 
sign  variations  in  the  first  column  of  the  Routh  Table  is  at  least 
N-l  (where  2(N-1)  is  the  highest  power  in  Y(iu)). 

If  all  of  the  necessity  tests  have  been  passed  but  no  suffi- 
ciency condition  has  been  satisfied,  a  modification  of  Siljak's  neces- 
sary and  sufficient  test  for  non-negativity  is  made  as  follows  [S3, 
K10]: 

1.  Number  the  rows  of  the  Routh  Array  starting  at  the  top 
r  =  1,2,. . . ,2N+1. 

2.  Denote  each  row  preceding  a  zero  row  by  r  ,V=l,2,...,m 

where  the  first  row  is  r   and  the  last  row   r 

o  m+1 

3.  V   is  the  number  of  sign  variations  between  two  consecutive 

zero  rows  r   and  r 

V      v-1 

4.  TT  =  2r   -  r    -r     -v+V        V  =  1  2     m 

V     v    V-1    v+1     V    V+l  '  

TT  ,  =  r    -  r  -v  ,. 

m+l    m+1    m    m+1 

2 

5.  If  the  sum  of  Tf  over  all  odd  u  is  zero,  then  Y(ou  )  is 

non-negative. 
Usage 

ITEST  =  NONNEG  (PSI , NDIM, IPRT) 

PSI  -  Vector  of  coefficients  of  real,  even  polynomial  to  be 
tested  (only  even  power  terms  are  included) . 

NDIM  -  Dimension  of  polynomial,  i.e. ,  N+l. 

IPRT  -  Flag  to  control  printing  of  Routh  Table 

(if  computed)  and  results  of  non-negativity  test. 
IPR"'  -   0   results  in  no  printing. 


Input 
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Output 

NONNEG  is  set  to  one  of  the  following  integer  values  depending 
on  the  results  of  the  non-negativity  test. 

0  -  Error  in  NONNEG,  no  conclusion  about  non-negativity  of  PSI . 

1  -  Constant  term  of  PSI  is  negative,  PSI  is  not  non-negative. 

2  -  Highest  power  term  of  PSI  which  is  non-zero  is  negative, 

PSI  is  not  non-negative. 

3  -  PSI  has  an  odd  number  of  sign  variations  between  consecu- 

tive coefficients  (Descartes1  rule  of  signs),  PSI  is  not 
non-negative. 

4  -  PSI  passes  all  of  the  necessity  tests  but  fails  Siljak's 

test,  PSI  is  not  non-negative. 

5  -  PSI  has  all  non-negative  coefficients,  PSI  is  positive 

definite. 

6  -  PSI  has  an  even  number  of  sign  variations  and  its  Routh 

Table  has  at  least  N-l  sign  variations  in  the  first  column, 
PSI  is  positive  definite. 

7  -  PSI  passes  Siljak's  test,  PSI  is  non-negative. 


Restrictions 

Calls  function  subprogram  NSIGNV. 


Identification  NSIGNV 

FUNCTION  NSIGNV  (X,B,L) 
Purpose 

To  determine  the  number  of  sign  variations  between  consecutive 
terms  of  a  real  vector. 
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Method 

The  L  elements  of  the  input  vector  X  are  examined  for  sign, 
a  corresponding  i  1  is  placed  into  working  vector  B  and  zero  elements 
of  X  are  disregarded.   The  L-M  (where  M  is  the  number  of  zero  elements 
of  X)  elements  of  B  are  sequentially  multiplied  and  every  B(I-l) ■ B(I)< 0 
is  recorded  as  a  sign  variation. 
Usage 

NSIGN  =  NSIGNV  (X,B,L) 
Input 

X  -  real  vector  to  be  tested,  at  least  L  elements  long 

B  -  real  working  vector,  at  least  L  elements  long 

L  -  number  of  terms  of  X  to  be  tested 
Output 

NSIGNV  is  set  to  the  total  (integer)  number  of  sign  variations 
between  L  consecutive  terms  of  the  vector  X. 

C. 4   Spectral  Factorization 

If  Y(to)  is  non-negative,  the  system  configuration  being  tested 
is  optimal  and  by  Theorem  3.2  a  performance  index  with  a  unity  rank 
positive  semidefinite  state  weighting  matrix  can  always  be  constructed. 
This  weighting  matrix  is  formed  from  the  outer  product  of  the  vector 
of  coefficients  of  the  Hurwitz  spectral  factor  of  Y(w)  ordered  with 
constant  term  first  (Section  3.5).   This  section  presents  an  efficient 
and  accurate  method  for  the  extraction  of  Hurwitz  spectral  factors  of 
non-negative  real  even  polynomials.   The  procedure  is  somewhat  unusual 
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in  that  no  roots  of  the  polynomials  are  explicitly  computed.   This 
procedure  is  implemented  by  subprogram  SPFCTR  and  supporting  subroutine 
QFACT. 

Identification  SPFCTR 

SUBROUTINE  SPFCTR  (PSI , FACTR , NDIM, IER) 
Purpose 

To  compute  the  Hurwitz  spectral  factor  of  a  real,  even,  non- 
negative  polynomial. 
Method 

A  factor  i|r(s)  of  the  input  polynomial  Y(U))  is  to  be  computed 

so  that 

N+1      2i 
Y(uj)  =  tCjco)  ^r(-jaj)  =   S   a.cu 

i=0   x 

and  i)r(s)  has  all  left  half-plane  zeros. 

2         2 

First,  the  formal  substitution  of  a   =  ou   into  Y(o)  )  reduces 

the  order  of  Y  from  2N  to  N  and  maps  the  zeros  of  \jr(s)  (and  conse- 
quently zeros  of  Y(u>))  from 

±  a   ±  j0,  to:  (a2  -  32)  ±  j2o'0  . 
An  approximation  P(a)  ,  to  a  quadratic  factor  of  Y(cr)  is  then 
generated  using  Bairstow' s  method  [H3] ,  Y  is  divided  by  P  and  its 
order  reduced  by  two.   P  becomes  a  factor  F  of  the  spectral  factor 
after  its  coefficients  have  been  altered  to  reflect  the  inverse  of 
the  mapping  described  above.   This  is  done  by  computing  the  discrim- 
inant of  P  and  by  examining  its  sign  along  with  the  signs  of  the 
coefficients  of  P  to  determine  the  necessary  transformation  to 
compute  the  coefficients  of  F. 
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The  steps  in  the  preceding  paragraph  are  repeated  and  the 
quadratic  factors  F  of  the  spectral  factor  are  multiplied  together 
to  produce  the  complete  spectral  factor.   This  method  avoids  the 
inaccuracies  involved  in  determining  the  2N  roots  of  Y(cu)  and  con- 
structing \jf(s)  from  the  left  half-plane  zeros. 
Usage 

CALL  SPFCTR  (PSI ,FACTR, NDIM, IER) 
Input 

PSI  -  Vector  of  coefficients  of  the  polynomial  which  is  to  be 

factored,  only  even  power  terms  are  stored,  i.e.  ,  nth  term 
of  the  vector  is  the  coefficient  of  cu       (double 
precision) . 
NDIM  -  Number  of  coefficients  of  PSI,  i.e. ,  N+l. 
Output 

FACTR  -  Vector  of  coefficients  of  the  Hurwitz  spectral  factor  of 
PSI  (if  non-negative)  stored  from  smallest  to  largest 
power  terms  (double  precision). 
IER  -  Integer  error  code: 

0  -  no  error 

1  -  maximum  number  of  attempts  at  finding  quadratic  factors 

of  Y(a)  exceeded  (presently,  2(n+1)  -  factorization 
aborted. 

2  -  order  of  input  polynomial  has  been  reduced  by  roundoff 

or  overflow  -  factorization  aborted. 

3  -  remaining  linear  factor  of  Y(a)  has  a  positive  zero,  no 

real  spectral  factor  can  be  generated  -  factorization 
aborted. 
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Restrictions 

If  the  input  polynomial  is  not  non-negative,  a  real  spectral 
factor  does  not  exist;  nonetheless,  the  algorithm  may  still  produce  an 
attempt  at  a  spectral  factor,  without  detecting  an  error.   The  results 
of  this  procedure  can  only  be  relied  upon  if  the  non-negativity  of  the 
input  polynomial  is  insured. 

Calls  subprogram  QFACT. 

Identification  QFACT 

SUBROUTINE  QFACT  (C, IC, Q.LIM, DNORM, IER) 
Purpose 

To  compute  and  divide  out  a  quadratic  factor  of  a  real  poly- 
nomial using  Bairstow' s  method. 
Method 

Let  the  input  polynomial  to  the  algorithm  be: 

N 

E 

i=0 


C(s)  =   E  c.s1 


2 

and  an  initial  guess  for  a  quadratic  factor  be:   s  -  vs  -  w.   By  divi- 
sion, a  quotient  and  a  remainder  are  computed, 

2 
C(s)  =  (s  -  vs  -  w)  B(s)  +  Rs  +  T, 

if   R  =  T  =  0,  the  initial  guess  was  correct  and  it  is  a  quadratic  factor. 
If   R^T^O,  then  a  sequence  of  v' s  and  w' s  is  generated  so  that  in 
the  limit   R  =  T  =  0.   Bairstow' s  method  is  simply  the  simultaneous  solu- 
tion of 

R(v,  w)  =  0 
and  T(v,  w)  =  0 

by  Newton1 s  method  [H3] . 
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Once  a  satisfactory  approximation  to  a  quadratic  factor  has 
been  obtained,  C  is  replaced  by  the  quotient  and  the  order  is  reduced 
by  two. 
Usage 

CALL  QFACT  (C, IC, P.LIM.DNORM, IER) 
Input 

C  -  vector  of  coefficients  of  the  polynomial  to  be  factored, 
stored  from  smallest  to  largest  power  term,  i.e. ,  nth  term 
of  vector  is  coefficient  of  s     (double  precision) 
IC  -  number  of  coefficients  of  C. 

P  -  vector  of  coefficients  of  an  initial  guess  for  a  quadratic 

2 
factor  (on  input):   s  +  P(2)s  +  P(l)  (double  precision). 

LIM  -  integer  specifying  the  maximum  number  of  iterations  allowed. 

Output 

C  -  vector  of  coefficients  of  the  polynomial  with  quadratic 

factor  removed,  the  last  coefficient  is  now  one  (double 

precision) . 

IC  -  number  of  coefficients  of  reduced  C,  i.e. ,  value  on  input 

less  two. 

2 

P  -  approximation  of  quadratic  factor:   s  +  P(2)s  +  P(l) 

(double  precision). 
DNORM  -  value  by  which  coefficients  of  C  on  input  have  been  normal- 
ized, if  other  than  one,  DXORM=C(IC)  (double  precision). 
IER  -  error  code  set  to  one  of  the  following  integer  values: 

0  -  no  error  encountered 

1  -  no  convergence  within  LIM  iterations 
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-1  -  input  polynomial  is  constant  or  undefined  -  or  an  over- 
flow occurred  during  normalization 
-2  -  input  polynomial  is  of  degree  one 

-3  -  no  further  refinement  of  the  approximation  to  a  quad- 
ratic factor  possible;  division  by  zero,  overflow,  or 
an  initial  guess  outside  of  the  region  of  convergence. 
Remarks 

Some  of  the  FORTRAN  code  was  taken  from  parts  of  subroutine 
DPQFB  of  the  IBM  System/360  Scientific  Subroutine  Package  [II, 
pp.  193-197]. 

C.5   Sample  Program 

OPTI  is  a  sample  main  program  which  demonstrates  how  a  given 
design  configuration  may  be  tested  for  optimal ity  with  respect  to  a 
quadratic  performance  index  and  how  appropriate  performance  indices 
are  computed  if  they  exist. 

The  program  accepts  as  inputs  the  coefficients  of  the  closed- 
loop  and  open-loop  characteristic  polynomials  of  the  system  under 
investigation.   The  stability  of  the  closed-loop  characteristic  poly- 
nomial is  not  tested.   It  is  assumed  that  the  closed-loop  character- 
istic polynomial  resulted  from  a  design  technique  which  insured 
stability;  if  the  closed-loop  system  is  not  stable,  the  results  of 
program  OPTI  are  not  supportable. 

The  even  polynomial , 

v(iu)  =  |cp  (juu)J   -  |cp(jcu)|  , 
is  formed  and  tested  by  subroutine  NONNEG  for  system  optimality. 
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If  the  system  is  not  optimal  for  a  positive  semidef inite  state  weight- 
ing matrix,  that  is,  Y(cu)  }t   0,  the  program  terminates  analysis  of  the 
present  system  and  interrogates  the  input  file  for  further  instructions. 

If  the  system  is  determined  to  be  optimal,  corresponding  posi- 
tive semidefinite  state  weighting  matrices  (Q)  of  the  performance 
index, 

00 


J=J(xQx+u)  dt, 


2 

are  computed.   If  Y(u>  )  has  all  non-negative  coefficients,  the  diagonal 

Q  is  printed.   In  either  case,  the  unity  rank  factor  of  Q,  that  is,  h 

T 
where   Q  =  hh  ,  is  computed  by  spectral  factorization  using  subroutine 

SPFCTR  and  printed. 

The  accuracy  of  the  factorization  is  checked  by  performing  the 
product , 

tt(uj)  =  t|;(ju))ijf(-juj)  , 
(where  ty  is  the  approximate  spectral  factor  of  Y)  and  determining  the 
average  magnitude  difference  in  the  coefficients  of  rr  and  Y.   The 
results  of  the  accuracy  test  are  also  printed. 

Several  tests  for  numerical  difficulties  are  made  in  the  course 
of  the  analysis  and,  if  a  problem  arises,  an  appropriate  explanation  is 
printed  and  the  input  is  queried  for  instructions. 
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Ir.put 


Each  input  card  image  must  be  in  the  following  format: 


Field 

Columns 

Type 

Use 

1 

1-3 

Alpha 

Instruction 

flag 

2 

4,5 

Integer 

Coefficient 

number 

3 

6-30 

Real 

Coefficient 

Stored  internally  in  double  precision  and  read  on 
format  D25. 16. 


Field  1  contains  (left  justified)  one  of  the  following 
instruction  flags: 

P  -  The  data  on  this  card  (and  all  following  cards  until 
another  flag  is  encountered)  describe   a  coefficient 
of  polynomial  cp  (open-loop  characteristic  polynomial). 
PK  -  Same  as  P  except  identifies  coefficients  of  cp   (closed- 
loop  characteristic  polynomial). 
*GO  -  Signals  the  end  of  data  for  a  given  system,  commands  the 

program  to  begin  analysis. 
END  -  All  program  execution  is  to  terminate. 

Field  2  contains  an  integer  (right  justified),  i,  indicating 
that  the  datum  in  field  3  is  the  coefficient  of  the  (i-1)  power  term 
of  the  polynomial  indicated  in  field  1  of  this  or  a  previous  card. 

Field  3  contains  a  coefficient  value  in  (either  single  or  double 
precision)  FORTRAN  real  number  format.   Only  non-zero  coefficients  need 
be  entered  (i.e.,  the  coefficients  are  zeroed  prior  to  reading  each  case) 
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Sample  Case 

The  output  of  a  sample  case  is  included  in  Section  C.6. 

C. 6   Subprogram  Listing 

In  order  to  facilitate  reference  to  this  program,  the  listing 
has  been  divided  into  functional  groups  of  subprograms.   These  func- 
tional groups  are  indexed  below. 

I.   Test  for  Non-negativity  of  Real,  Even  Polynomials 

1 .  NONNEG 

2.  NSIGNV 

II.   Spectral  Factorization 

1 .  SPFCTR 

2.  QFACT 

III.   Computation  of  Magnitude  -  Square  Polynomials 
1.   MAGSQ 
IV.   Sample  Program 

1.  OPTI 

2.  Sample  Problem 
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FUNCTION  NONN£G( A2K.NDIMA,  IPRT)  NONN   ID 

C      FUNCTION  NONNEG  DETERMINES  THE  NON-NEGATIVITY  OF  EVEN  REAL  POLYS.  NUNN   23 

C  NUNN   30 

C         NONNEG  RETURNS  THE  FOLLOWING  COOESO  NONN   40 

C         0  -  ERROR  IN  NONNEG  NUNN   50 

C          1  -  CONSTANT  TERM  NEGATIVE  NUNN   60 

C          2  -  HIGHEST  POWER  (NONZERU)  TERM  NEGATIVE  NUNN   70 

C          3  -  POLYNOMIAL  HAS  ODD  NUMBER  OF  SIGN  VARIATIONS  NUNN   30 

C          4  -  POLYNOMIAL  HAS  EVEN  NO.  OF  SIGN  VARS.  HAS  LESS  THAN  N-l  NUNN   90 

C              SIGN  VARS.  IN  1ST.  COL.  OF  ROUTH  TABLE  FAILS  SILJAK'S  TEST  NUNN  100 

C         5  -  POLYNOMIAL  COEFICIEnTS  ALL  POSITIVE  NONN  110 

C          6  -  POLYNOMIAL  HAS  EVEN  NU.  OF  S IGN  VARS.  ROUTH  TABLE  HAS  AT  NONN  120 

C             LEAST  N-l  SIGN  VARIATIONS  IN  1ST.  CULUMN  NONN  130 

C          7  -  POLYNOMIAL  HAS  EVEN  NO.  OF  SIGN  VARS.  RUUTH  TABLE  DUES  NOT  NONN  140 

C              HAVE  N-l  SIGN  VARS.  IN  1ST.  COLUMN  PASSES  SILJAK'S  TEST  NONN  15J 

C                                                             2 (N-l I  NONN  163 

C      A2K  -  INPUT  POLYNOMIAL  COEFS.   A2<(N)  =  COEF.  OF  W  NONN  1?0 

C      NDIMA  -  NUMBER  OF  NONZERO  COEFICIENTS  NONN  l-lO 

C      IPRT  -  FLAG  10  INDICATE  WHETHER  THE  NORMALIZED  ROUTH  TABLE  AND  NONN  193 

C              THE  RESULTS  OF  NONNEG  ARE  TO  BE  PRINTED  (IPRT.Nc.O)  NONN  200 

C                                                  J.M.  ELDER   4/23/71  NONN  210 

DOUBLE  PRECISION  DR1  ,.DR2,  DTEMP  ,  R  1M  ,  R2M  ,  SIGNR  ,RTSC  NONN  220 

DIMENSION  A2K(  1  )  ,R( 40, 20),R1(  1  ) ,0<  1 ) »  JR0WI40) ,NV( 1 )  ,  NONN  230 

IDRH20),DR2(20)  ,DTEMP(  20),  IUUT(  12,  8)                      .  NONN  243 

COMMON  R.DR1.DR2, DTEMP, RIM, R2M  NONN  253 

EQUIVALENCE  ( R (  1 ,1  )  , R  1  (  1 )  ) ,  ( R (  1 , 2 ) , NV(  1  ) ) ,  ( R ( 1 , 3 ) , 3 (  1 )  )  NONN  260 

DATA  I0UT/4H.N0T.4H-J0  C,  4HUNCL  ,  4HUS  10,  4HN  PO  ,  4HSS  I  3  i4H  LE  ,  ,  4HE  RRO  ,  NONN  270 

14HR  INf4H  NUN.4HNEG  , 4H     , 4H . NOT , 4HC0NS , 4H TAN T , 4H  TER.4HM  OF,  NUNN  290 

24H  PUL,4HYN0M,4HIAL  , 4HNEG A , 4HT I VE , 2»4H     , 4H . NOT , 4HH I GH , 4HE i T  ,  NONN  290 

34HP0WE.4HR  C0,4HEFFI,4HCItN, 4IIT  NE , 4HGAT I , 4H VE   ,2»4H     .4H.NOT,  NUNN  300 

44HP0LY,4HN0MI,4HAL  H.4HAS  0,4r,DD  N.4HJ.  0,4HF  SI.4HGN  V.4HARIA,  NONN  310 

54HTI0N.4HS    ,4H.NUT,4HPASS,-4HES  ^ , 4HEC E S , 4H SI T Y , 4H ,  3U.4HT  FA,  NONN  320 

64HILS  ,4HSJFF,4HICIE,4HNCY  , 4HT EST , 4H. . . . , 4HP3L Y , 4HN0M I , 4H AL  H,  NONN  330 

74HAS  A.4HLL  P, 4 HOS I T , 4H I VE  , 4HC0EF , 4HF I C I , 4HEN T S , 4H     ,4H NONN  340 

84HR0UT,4HH  TA.4H8LE  ,4HHAS  ,.4H  N  S.4HIGN  .4HVARS.4H.  IN.4HFIRS,  NUNN  350 

94HT  CU,4HL.   ,  4H.  .  .  .  ,.4HP0L  Y  ,  4HNOM  I  ,  4HAL  P.4HASSE.4HS  SU,4HrFIC,  NONN  360 

A4HIENC.4HY  TE.4HST   ,2«4H     /  NONN  370 

DATA  IMIN/1H-/, IZER/1H0/, IPLUS/1H*/  NUNN  380 

10  FORMAT  (///,  2GX  ,  22HNORM  AL  I  Z  ED  RUUTH  ARRAY/6X.6H  ROW   ,  5  X ,.  NONN  390 

115HNORMALIZED  RUWS/1H  , 70I1H-))  NONN  400 

20  FORMAT (6X.1H  ,I3,2H   ,  1 Al ,  1P8E  1 1 . 2 )  NONN  410 

33  F0RMATI6X.1H  ,I3,2H   ,  1  A  1 ,  1  P8E  1 1 . 2 , t /6X , 1H  ,4X,2H   .1P8E11.2))  NUNN  420 

40  FORMATI3H  »R,I2,2H   ,I3,2H   ,  1 A  1  ,  1  P8E  1 1  .  2  )  NU.'I.M  430 

50  FURMAM3H  »R,I2,2H   ,I3,2H   ,  1 A I ,  IP  8F  1 1  .  2  ,  NONN  440 

K/6X.1H  ,4X,2H   .1P3EU.2))  NUNN  450 

60  FORMATI28H  »  ROWS  PRECEEDING  ZERO  R0WS/66H   (ZERO  ROWS  HAVE  BEEN  RNONN  460 

1ESOLVED  BY  DIFFERENTIATING  PRECEEDING  ROW))  NONN  470 

70  FORMAM  49X.1H2/32H   PJLYNUMIAL  HAS  BEEN  MULTIPLIED,  NONN  483 

119H  BY  THE  FACTOR   S  +.F7.4)  NONN  493 

80  FORMAT(/10XflH2/5X,12HPSI( W  )  1S...1A4.18H  POSITIVE  SEMI  DEF,  NONN  500 

1  7HINITE,  ,  11A4)  NUNN  510 

ICHNG=0  NONN  523 

C      A2M 1) ...CUNSTANT  TERM  MUST  BE  NON-NEGATIVE  FOR  NON-NEGATIVITY  NONN  533 

IF( A2K(  1)  )90,10O, 100  NUNN  540 

90    N0NNEG=1  HQXU    550 

GU    TO       630  NUNN    560 

C               REDUCE    DIMENSION    IF    HIGHEST    ORDER    TERMS    NEGLIGIBLE  NUNN    570 

100     1F(ABS( A2KINDIMA) ) -1 .0 E-50 ) 1 10 , 110,120  NUNN    580 

110    NDIMA=NDIMA-1  NCNN    59° 

GO    TO       100  NOfJN    600 

C               TEST    FOR    HIGHEST    ORDER    TERM    POSITIVE  NONN    610 

120    NDIK*NOIMA  NUNN    620 

IF ( A2K( NO  I M) 1130,6  20, 140  NUNN    6  30 

130    N0NNEG=2  HOfiU    6''3 


135 


GO  TO   630  NONN  653 

C      CHECK  FOR  POSUIVITY  USING  DESCARTES'  RULE  OF  SIGNS  NUNN  653 

C          DETERMINE  NUMBER  OF  SIGN  VARIATIONS  NONN  670 

1*0  ICHNG=NSIG-IV(A2K,B,NDIM)  NONN  683 

IF( ICHNG)620,150, 160  NONN  690 

150  NONNCG=5  NONN  703 

GO  TO   630                   ■  NOUN  710 

160  FICHNG=FLOAT(lC  '  G)/2.0*1.0E-*  NONN  723 

C      TEST  FOR  ODD  ML   ER  OF  SIGN  VARIATIONS  NONN  733 

IF(2»INT(F1CHNG/-ICHNG) 170,  100,620  NONN  7*0 

170  N0NNEG=3  NONN  750 

GO  TO   630  NONN  763 

C      COMPUTE  ROUTH  ARRAY  NONN  773 

180  DO   190  I=1,NDIM  NONN  780 

190  DTEMPI I )=DBLE< A2K( I ) )  NONN  793 

C      FILL  FIRST  TWO  ROWS  NONN  803 

NCOL10=0  NONN  313 

200  IZER0=1  NONN  820 

NROWI 1 )=1  NONN  830 

S1GNR=  1.000  NONN  BO 

J=NDIM*1  NONN  853 

R1M=O.ODO  NONN  860 

R2M=.O.OD0  NONN  870 

DO   2*0  I=1,NDIM  NONN  e83 

J»J-1  NONN  690 

0R1(J)=SIGNR»DTEMP( I)                                     .  NONN  900 

DR2IJ)=2.030*DBLE( FLOAT! 1-1) )»DR1( J)  NONN  910 

IFlDABStDRl (J) J-R1M1220, 220,210  NONN  923 

210  R1M=DABS(D31 I  J)  )  NONN  933 

223  IF(  DAbStDR2  (J)  )-R2M)2*0,  2*0,230  NOvIN  90 

230  R2M=DABS<  DR2U  )  )  NONN  953 

2*0  SIGNR=-SIGNR  NONN  963 

C      NORMALIZE  THE  FIRST  TWO  ROWS  NOUN  973 

DO   250  J=1,NDIM  NUNN  9C0 

0R1 ( J)=0R1( JJ/R1M  NONN  990 

DR2( J)=0R2( J >/*2M  NONN1C00 

R(1,J)=SNGL!DR1( J) )  NONN1010 

250    RI2, J)=SNGL(DR2! J)  >  NUNN1020 

C               COMPUTE    REST    OF    ARRAY  NOUN1030 

NEND1  =  2»NDIM-1  NONNK*0 

DO      *00    l=3*NEN01  N0lNij53 

NEND  =  NDIM-l NT (<  FLOAT!  I  1-1.0 )/2. 0  +  1. OE-*)  NUNN1060 

SIGNR=1.0DO  NUNN1070 

IF(DR2(  1 J  1260,260, 270  N0NN1CS3 

260  SIGNR=-1.0D0  NONN109J 

270  R1M=0.0U0  NQNN1100 

DO   290  J=1,NEND  NUNN1110 

DTEMPl  J)=SIGNR«(0R2(  ll»OMI  J  +  IJ-DU!  1>»DR2!  J  +  l)  )  NONN1123 

IF (DABS (DTE  MP ( J) ) -RIM ) 290, 2  90, 230  N0NN1133 

280  R1M=DABS(DTEMP(J) )  NONN11*0 

290  CONTINUE  NONN1150 

C      CHECK  FOR  ZERO  ROWS  IN  ARRAY  NONN1160 

IFIR1M  -1.30-12)300,300,3*0  NONN1173 

C      NOTE  THAT  ZERO  ROW  HAS  BEEN  ENCOUNTERED  AND  STORE  ITS  ROW  NO.  -1   NUNN i  183 

300  IZERO=IZERO+1  NUNN1 1 90 

NROWI IZfcRO) =1-1  NONNl?r0 

C      RESOLVE  ZERO  ROWS  N0NM213 

NPWR»2.NDIM-I  NONN1223 

ftlM-O.ODO  NUNN1233 

00   330  J=1,NDIM  NONN12*0 

IF(NPWR)3*0,.3*0,31O  NUNN1253 

310  DTEMPt J)xOBLE(FLOAT(NPWR)).DR2(  J)  NUNN1260 

IF(DABSIDTEMP(  J  )  )-RlM)330,  330,320  NUWN12  70 


320  R1M=DABS(DTEMP( J) ) 
330  NPwR=NPwR-2 

NORMALIZE  ROW 


NONNi; 3  3 
NONN 12 93 
NUNN1300 
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3*0    00       350    J=l,.NEND  N0NN1310 

DR1 ( J)=DR2< J )  NU.NN132D 

0R2UJ=DTEMP(  J)/RIM  NONN1330 

350    R(  I  ,J)=SNGL(DR2U)  )  NUNN1H0 

C      TEST  FOR  ZERO  IN  FIRST  COLUMN  NUNN1350 

IF<DABS!DR2(  1)  )-l  .00-12  )  360,  3.60,  390  NONN1360 

C      MULTIPLY  IN  AN  OFFSETTING  FACTOR  AND  BEGIN  TA3LE  OVER  NONN1370 

360  IF(NC0L10-10)370,620,620  NUNN1380 

370  RTSC=1.  11  lllllllll  1111D0+D3LE!  FLOAT!  NCOL  10  )  )  NONN1390 

NDIM=NDIMA*1  NONN1400 

DTEMP! 1  )*08LE(A2K(  1)  )*RTSQ  N0NN1413 

DTEMP (NDI Ml =D3LE<  A  2K I  NO  IMA)  )  N0NN14  20 

00   380  I  =  2».NDIMA  N0NN1430 

380  DTEMP! I )=RTSQ»DBLE(A2K( I ) )+DBLE(A2K( 1-1) 1  NUNN144D 

NCOL10  =  .NCOL10*1  NONN1450 

GO  TO   200  NUNN1460 

390  NC0L10  =  0  NUNN1473 

DR1 ( J*i )=DR2( J+l)  NUNN1483 

DR2 ( J+1)=O.ODO  NUNN1490 

400  RU,JH)»0.0  NONN1500 

C      INCLUDE  LAST  ROW  N0NN1510 

IZERO=  IZcRO+.l  NUNN1520 

NROW! IZER0)=NEND1  N0NN153D 

IF(  IPRT)413,.530,410  NUNN1540 

C      PRINT  NORMALIZED  RUUTH  TABLE  NU.NN1550 

410  WRITEI6,   10)  NONN1563 

KLESS1=0                                  '  NONN1570 

K=l  NUNN1580 

DO   520  I  =  l,.NENDl  NUNM1590 

IF(R(  1,1)  1420,430,440  NUN,<I1600 

42  0  ISIGN=IMIN  N0NN161D 

GO  TO   450  NUNN1623 

430  ISIGN=IZER  NUNN1630 

GO  TO   450  NO.NN1640 

440  ISIGN=IPLUS  NONM653 

450  NENO=NDIM-INT(  (  FLOAT!  I  1-1.01/2. 0+:i.0E-4)  N0NN1663 

IF  (  I -NROW (K  1)490,460,490  N0.NN16  70 

460  K-K+l  NUriN1680 

KLESS1=KLESS1+1  N0NN1693 

IF!NEND-8)470,470,480  NONN1703 

470  WRITEt6f   401KLESS1,  I,  ISIGN,.(R(  I..J  1  ,J=1,NEND)  N0*N1710 

GO  TO   520  NONN1720 

480  WRITE!6,.   50 )KL ESS  1, I ,  I S IGN , ( R ( 1 , J ) , J  =  1 , NEND )  NJNN1733 

GO  TO   520  MONN1740 

490  IFtNEND-81500,500,  510  N0NN1750 

500  WRITE16,.   20)1, IS1GN, (R! I, J),J=1,NEN0)  NOUNl7bO 

GO  TO   520  NO.NN1770 

510  WRITE16,.   30)1,  ISIGN,  (R(  I,  J  >  ,.J  =  1 ,  NEND)  N0MN1780 

520  CONTINUE  N0NN1790 

FICHNG=SNGL!RTSO)  NUNN1800 

IF(NDIM.GT.NDIMA)WRITE!6,2010)FICHNG  NUNN1813 

WRITE16,   60)  NON.N1820 

C      DETERMINE  NUMBER  OF  SIGN  VARIATIONS  IN  FIRST  COLUMN  NUWN1830 

530  ICHNG=NSIGjV!R1,B,NcND1 )  NUNN1840 

IF! ICHNG-N0IM-2)550,540, 540  NONN18  50 

C      SIGN  VARIATIONS  SHOULD  BE  AT  LEAST  N-l  FOR  POSITIVITY  N0MN1863 

540  N0NNEG=6  NUNN1870 

GO  TO   630  N0N.U800 

C  N0NN1H93 

C      SILJAK'S  TEST  FOR  NON-NEGATIVITY  NUNN1900 

C  NUNN1910 

550  00   560  t  =  2,-IZER0  NUNN1920 

J=NROWl I J-NROW! 1-1 )  N0NN1933 

K=NROW( 1-1) »1  NUNN1940 

560  NV(  I )=NSIGNV(Rl(K),B,JI  NUN. 1950 

C      TEST  UDO  SECTOR  NUMBERS  NUNN1963 
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NPI=0  NUNN1970 

DO   570  I=1,IZER0,2  N0NN1963 

IF(  I  +  l-IZERO)570,580,.593  N0NN1993 

5  70  NPI=2«NR0W(  I  ♦  1  )-N«t)W(  I*2)-NR0W(  I  )  ♦  2»  {  N  V  t  1  ♦  2  )  -N  V  (  I  ♦  1  1  )  «•  NPI  NUNN20O3 

580  NPI=NROW( UERd)-NROW( IZER0-l)-2'NVl IZE201+NPI  N0NN2G13 

590  IFtNPI  )610,600t610  NONN2020 

600  N0NNEG=7                        '                              .  NONN2033 

GO  TO   630  NUNN2040 

610  N0NNEG=4  NONN2050 

GO  TO   630  NONN2063 

620  NONNEG=0  N0NN2070 

630  IF(  1PRT)643,650,640  NONN2OS0 

640  J=.N0NNEG*1  N0NN2093 

WRITE16,.   80)I0UTl  1,  J  1  .  I  IOUT(  I, J),  1  =  2,  12)  NUNN2100 

650  RETURN  NONN2110 

END  NUNN2120 


FUNCTION  NSIGNVIX.W, LENGTH) 
C      FUNCTION  NSIGNV  COMPUTES  THE  NUMBER  OF  SIGN  VARIATIONS 
C  IN  A  VECTOR  Xt  LENGTH  LONG 

C  W  IS  A  SINGLE  PRECISION  WORK  VECTOR  OF  LENGTH    LENGTH 

C  NSIGNV  IS  USED  IN  CONJUNCTION  WITH  FUNCTION  NOMNEG 

DIMENSION  XI  1),W(  1  ) 
C      COUNT  NUMBER  OF  SIGN  VARIATIONS 

NSIGNV=0 

NDIM=0 

DO    20  1  =  1, .LENGTH 

If  U<  I  )  )  10,20,  10 
10  NDIM=NDIM+1 

W(NDIM)=SIGNI1.0,X(  I)  ) 
20  CONTINUE 

DO    40  I=2,NDIM 

IF(W( I-1)»W( I) )30, 40,40 
30  NSIGNV=NSIGNV+1 
40  CONTINUE 

RETURN 

END 


NSIN 

10 

NSIN 

20 

NSIN 

3D 

NSIN 

40 

NSIN 

50 

NSIN 

60 

NSIN 

73 

NSIN 

faO 

NS  I  N 

9J 

NSIN 

100 

NSIN 

110 

NSIN 

123 

NSIN 

130 

NSIN 

140 

NSIN 

153 

NSIN 

160 

NSIN 

170 

NSIN 

180 

NSIN 

190 

NSIN 

200 
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SUBROUTINE  SPFCTR!  PS  I , F ACTR , ND  IMA,  IER) 

SU3ROUTINE  SPFCTR  COMPUTES  THE  SPECTRAL  FACTOR  X(JW)  OF 

2 
AN  EVEN  POLYNOMIAL  Z(W  ),I.E. 


2 

2(W  ) 


XI  JW) 


=  X( JW)X( JW) 


A2K 


2  2IN-1) 

INPUT  POLYNOMIALS!  W  ))   A2K{N)  =  COEF.  OF  W 


U 
20 
30 
t*0 
53 
60 
70 
80 
93 


N-l 


RE 


20 


FACTR=  SPECTRAL  FACTOR  POL YNOM I AL ( X t S ) )   FACTR(N)  =  COEF.  OF  S 
NDIMA=  NUMBER  OF  COEFICIENTS  OF  FACTR  (ORDER  ♦  1) 
IER   =  ERROR  CODE 

IER  =  0   NO  ERROR 

IER  =  1   MAX.  NO.  OF  ATTEMPTS  AT  FINDING  QUAD.  FACTORS 

EXCEEDED 
IER  =  2   ORDER  OF  INPUT  POLYNOMIAL  HAS  BEEN  REDUCED  BY 

ROUNDOFF  OR  OVERFLOW  OCCURRED  DURING  NORMALIZATION 
IER  =  3   REMAINING  LINEAR  FACTOR  HAS  POSITIVE  ZERO 

MARKSO 

1.  IF  FACTORIZATION  IS  SUCCESSFUL ( I ER  =  0),  FACTR  MAY  STILL  NUT 
BE  THE  SPECTRAL  FACTOR  IF  THE  INPUT  POLYNUMIAL  IS  NUT  NON- 
NEGATIVE  DEFINITE   IT  HILL  DIFFER  FROM  THE  SPECTRAL  FACTOR 
IN  THAT  THE  RIGHT  HALF  PLANE  ZEROS  OF  THE  SPEC.  FACTOR  WILL 

BE  REFLECTED  ABOUT  THE  IMAGINARY  AXIS.   FACTR  WILL  ALWAYS  HAVE 
LEFT  HALF  PLANE  ZEROS. 

2.  SPFCTR  USES  SUBROUTINE  U'FACT  TO  DETERMINE  QUADRATIC  FACTORS. 

3.  PSI  AND  FACTR  ARE  DOUBLE  PRECISION  ARRAYS 

J.M.  ELDER   4/23/71 
DOUBLE  PRECISION  F ACTR ( 1 ) , PS  I ( 1 ) , WORK! 20) , QUAD { 4) , 
1A0,A1,R0,R1,F0,F1,DISC,DNDRM 
EQUIVALENCE  I AO , QU AD(  1 )  ) , ( A 1, QUAD! 2) ),  I RO, QUAD (3) )  , (Rl  ,QUAD(4)  ) 
COMMON  W0RK,QUAD,F0,F1,DISC,STRT,DEI 
NDIM=NDIMA 
DNORM=1.0D0 
DO    10  I  =  UNDIM 
FACTRI  I  )=0.0D0 
SET  NMAX.MAX.  NO.  OF  TRYS  AT  DETERMINING  ALL  THE  FACTORS  OF  PSI 

(PRESENTLY  TWICE  NDIM)  AND  I TM*X , MAX. NO.  OF  ITERATIONS  PER  T^Y 
NMAX=2»NDIM 
NTRY=1 
ITMAX=100 
FACTRt 1 )=1.0D0 
NDIMF=l 
STRT=1.0E*4 
DEL=STRT/FLOAT(NMAX) 
GENERATE  A  NEW  GUESS..IF  REQUIRED  (Al  0,T0  TRY  TO  FORCE  REAL 

DISTINCT  ZEROS  OF  QUADRATIC  FACTOR  OF  COMPANION  POLYNOMIAL) 
A0=DBLE(STRT) 
A1=2.5D0»DBLE(STRT) 
STRT=STRT-DEL 

HAS  MAXIMUM  NUMBER  OF  ATTEMPTS  AT  FACTORING  PSI  BEEN  EXCEEDED 
IF (NTRY-NMAX 14  0,40,340 
NTRY=NTRY+1 

APPROX.  QUAD.  FACTOR  OF  COMPANION  POLY.  (PSI)  USING  BAIRSTUWS  M. 
CALL  QFACTl PSI ,NDIM,QUAD, I TMAX,DN3RM, IERQ) 
TEST  FOR  ERRORS  IN  COMPUTING  FACTOR 
IF( IERQ)50,60,2O 
IF(  IER0«-.2)20,26O,.350 

MULTIPLY  FACTR  BY  LHP  QUADRATIC  PART  OF  QUARTIC  FACTOR 
COMPUTE  LHP  QUADRATIC  PART  COEFFICIENTS 


SFTR 
SFTR 
SFTR 
SFTR 
SFTR 
SFTR 
SFTR 
SFTR 
SFTR 
SFTR  100 
SFTR  110 
SFTR  120 
SFTR  133 
SFTR  143 
SFTR  150 
SFTR  160 
SFTR  170 
SFTR  180 
SFTR  lf,0 
SFTR  200 
SFTR  210 
SFTR  220 
SFTR  230 
SFTR  240 
SFTR  250 
SFTR  260 
SFTR  270 
SFTP  280 
SFTR  290 
SFTR  303 
SFTR  310 
SFTR  320 
SFTR  330 
SFTR  340 
SFTR  350 
SFTR  360 
SFTK  370 
SFTR  380 
SFTR  393 
SFTR  400 
SFTR  413 
SFTR  420 
SFTR  430 
SFTR  *W 
SFTK  450 
SFTR  460 
SFTR  470 
SFTR  480 
SFTR  493 
SFTR  500 
SFTR  510 
SFTR  520 
SFTR  533 
SFTR  540 
SFTK  550 
SFTK  560 
SFTR  570 
SFTR  580 
SFTR  590 
SFTR  600 
SFTR  610 
SFTR  620 
SFT?  630 
SFTP  640 
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„  •  SFTR  650 

60  F1=0.0D0  SFTR  66D 

C      TEST  DISCRIMINANT  TO  DETERMINE  WHERE  ROOTS  FALL 
DISC= Al«»2-4.0D0»A0 
IF ( DABS( 0  ISC 1-l.OD- 10 ) 70,70,80 

70  oisoo.ooo  "':H'  ;~: 

80  1F(DISC>190,170,9Q  ■  *Ll"  «*" 

C     PAIR  OF  REAL  ZEROS 
90  IF (AO) 100, 130, 130 
C     ZEROS  IN  OPPOSITE  HALF  PLANES 

C         MULTIPLY  IN  EXTRA  FACTOR  »£IJJ  ,JJ 

100  AO=DSURT(DISC) 

F0=0.5  00«DABS(A1+A0I 
JMP»0 

GO  TO   280 

C         DIVIDE  REDUNDANT  FACTOR  OUT  OF  PSI  ^MK  dOJ 

HO  F0  =  -0.5D0»(A1-A0) 
RO=PSl  (NDIM). 
NDIM=NDIM-1 
J=NDIM 
120  R1=PSI(J) 
PSI (J)=R0 
RO=R1-RO«FO 

IF(J)200.200,120  |FTR  890 

C      REAL  DISTINCT  ZEROS  IN  SAME  HALF  PLANE  •          SFTR  900 

130  IFIAIU^O.190,190  *FTR  910 

C      REAL  DISTINCT  ZEROS  IN  RIGHT  HALF  PLA-IE  SFTR  920 

1*0  IF(JMP-2)150,160,150  SFTR  930 

150  AO  =  DSURT(UISC)  »"R  9*0 

F0=-0.5D0»l Al-AO)  bM  K  V5J 

JMP  =  2  SFTR  9(>J 

GO  TO   210  l™    l10 

160  F0  =  -0.5D0MAl  +  A0).  SMK  VUJ 

GO  TO   200  •                    5f"TK  ^'J 

C      ZERO  UF  MULTIPLICITY  TWO  SFTR1000 

170  F0=0.5D0»DABS(A1>  lelolnon 

IF(A1)210, 2n, 180  ccrolr^ 

180  F1  =  DSQRT(2.0D0»A1)  SUKUJJ 

GO  TO   200  SFTRIj^O 

C      PAIR  OF  COMPLEX  ZEROS  OR  REAL  ZEROS  IN  LEFT  HALF  PLANE             SFTR1050 

190  FO=DSORT(A0)  SFTRI060 

F1=DSQRT(A1+2.0DO»FO)  SFTRU70 

C      MULTIPLY  IN  NEW  SPECTRAL  FACTOR                                     SFTR1P30 


SFTR  680 
SFTR  690 
SFTR  700 


SFTR  730 
SFTR  740 
SFTR  750 


SFTR  770 
SFTR  780 
SFTR  790 


SFTR  820 
SFTR  830 
SFTR  8'.0 
SFTR  850 
SFTR  660 
SFTR  870 
SFTR  880 


200  JMP=1  SFTR1090 

210  NDIMF=NDIMF+2  SFTR1100 

FACTR(NDIMF)=1.0D0  SFTR1H3 

NEND*NDIMF-3  SFTR1120 

IF(NEN0)2*0,.2*0,.220  SFTR1130 

220  DO   230  I=1,NEND  SFTR1140 

J=NDIMF-I  SFTR1150 

230  FACTRt J)  =  FACTR( J-2 ) +F 1 • F ACTRI  J-l )*F0*F ACTR  (  J  )  SFTR1160 

2*0  FACTR(2)=F0»FACTR(  2)+Fl»FACTR(  1)  SFTR1170 

FACTRl 1)=FD»FACTR( 1)  SFTRUBO 

IS  SPtCTRAL  FACTOR  COMPLETE  SFTR1190 

IFINDIM-1 )250, 320, 250  SFTR1200 

2'   GO  TO  (20,30),JMP  SFTR1210 

MULTIPLY  IN  FIRST  ORDER  FACTOR  SFTR1220 

260  FO=PSI  (  1  )  SFTP.123D 

270  IF(F0)360,23O,28O  SFTR1240 

280  FO  =  DS(jRT(F0)  SFTR125J 

NDIMF  =  NOIMF+:i  SFTR1260 

FACTR«NDIMF)»=1.0D0  SFTR1270 

NEND=NDIMF-2  SFTR12B0 

1F(NEND)310,.310,290  SFTR1290 

290  DO   300  I  =  UNEND  SFTR1300 
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J=NDIMF-I  SFTR1310 

330  FACTR1 J)=F0»FACTR( J  J f FACTK ( J-l )  SFTR1320 

310  FACTRf 1)=F3«FACTR( II  SFTR1330 

IF(JM|')320,  110,320  SFTR13A0 

IER=0  SFTR1353 

UNNORMALIZE  FACTR  AND  PSI      •  SFTR1360 

DISC=DSQRTIDNORM)  SFTR1370 

DO   330  I=1,NDIMA  SFTR13S0 

FACTRU)  =  DISC»FACTR(  I  )  SFTR1390 

PSI ( I )=DNORM»PSI(l)  SFTRlsOO 

RETURN  SFTR1410 

MAXIMUM  NUMBER  OF  ATTEMPTS  AT  FACTORING  PSI  EXCEEDED  SFTR1<»20 

IER=1  SFTRU30 

RETURN  SFTR1<,<,0 

ORDER  OF  LAST  FACTOR  OF  PSI  REDUCED  BY  ROUNDOFF  OR  OVERFLOW  SFTRU50 

OCCURRED  DURING  NORMALIZATION  OF  POLYNOMIAL                     SFTR1460 

IER=2  SFTR1470 

RETURN  SFTR1480 

REMAINING  LINEAR  FACTOR  OF  COMPANION  POLY.  HAS  POSITIVE  ZERO  SFTRl^.90 

IER=3  SFTRibOO 

RETURN  SFTR1513 

EMO  SFTR1520 


320 


330 


3^0 


350 


360 
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SUBROUTINE  QFACTIC, IC, Q , L IM , DNORM, IER)  OFAC  13 

C      SUBROUTINE  QFACT  COMPUTES  AND  REMOVES  OUADRATIC  FACTORS  FROM  QFAC  20 

C          REAL  POLYNOMIALS  USING  BAIRSTOW'S  METHOD.    IT  IS  A  MODI F I C AT  I ONCF AC  30 

C         OF  SUBROUTINE  DPOFB  OF  THE.  SYSTEMS  360  SCIENTIFIC  SUBROUTINE  QFAC  40 

C         PACKAGE.                                                •  0FAC  53 

C             C      DOUBLE  PRECISION  INPUT  VECTOR  CONTAINING  THE  OFAC  60 

C                   COEFFICIENTS  OF  P(X)  -  C(l)  IS  THE  CONSTANT  TERM  QFAC  70 

C                   (DIMENSION  IC)                                  -  CFAC  80 

C             IC   -  DIMENSION  OF  C  QFAC  90 

C             0      DOUBLE  PRECISION  VECTOR  OF  DIMENSION  4  -  ON  INPUT  OIDCFAC  100 

C                   AND  0(2)  CONTAIN  INITIAL  GUESSES  FOR  01  AND  U2  -  ON  CFAC  110 

C                   RETURN  0(1)  AND  0(2)  CONTAIN  THE  REFINED  COEFFICIENTS  CFAC  120 

C                    Ql  AND  02  OF  Q(X),.  WHILE  0(3)  AND  0  ( 4 )  CONTAIN  THE  OFAC  130 

C                   COEFFICIENTS  A  AND  B  OF  A+B»X,  WHICH  IS  THE  REMAINOER  OFAC  140 

C                   OF  THE  QUOTIENT  OF  P(X)  BY  Q(X)  QFAC  150 

C             L1M  -  INPUT  VALUE  SPECIFYING  THE  MAXIMUM  NUMBER  OF  QFAC  160 

C                   ITERATIONS  TO  BE  PERFORMED  OFAC  170 

C             DNORM-NONZERO  COEFFICIENT  OF  HIGHEST  POWER  TERM, BY  WHICH  CFAC  lbO 

C                   C  IS  NORMALIZED  QFAC  190 

C             IER  -  RESULTING  ERROR  PARAMETCR  (SEE  REMARKS)  OFAC  200 

C                   IER=  0  -  NO  ERROR  OFAC  210 

C                   IER=  1  -  NO  CONVERGENCE  WITHIN  LIM  ITERATIONS  QFAC  220 

C                   IER=-1  -  THE  POLYNOMIAL  P(X)  IS  CONSTANT  OR  UNDEFINED  QFAC  230 

C                             -  OR  OVERFLOW  OCCURRED  IN  NORMALIZING  PU)  CFAC  240 

C                   IER=-2  -  THE  POLYNOMIAL  P(X)  IS  OF  DE-GREE  1  QFAC  250 

C                   IERt-3  -  NO  FURTHER  REFINEMENT  OF  THE  APPROXI MATI UN  TOCFAC  260 

C                             A  OUADRATIC  FACTOR  IS  FEASIBLE,  DUE  TO  E1THER0FAC  270 

C                             DIVISION  BY  0,  OVERFLOW  OR  AN  INITIAL  GUESS  OFAC  280 

C                             THAT  IS  NOT  SUFFICIENTLY  CLOSE  TO  A  FACTuR  OFCFAC  290 

C                           P<X)  QFAC  300 

C                                              J.H.  ELDER  4/29/71  QFAC  310 

DIMENSION  C(  1), 0(1)  QFAC  320 

OOUBLE  PRECISION  A, 3, AA, B8 , CA, CB.CC ,CD , Al , 3 1 ,C 1 ,H,HH ,01 , 02 , 001  ,  OFAC  330 

1                  QQ2,.Q0Q1,CQQ2,DQ1,DQ2,EPS,EPS1,C,Q,DNURM  QFAC  340 

C  QFAC  350 

C         TEST  ON  LEADING  ZERO  COEFFICIENTS  QFAC  360 

IER^O  OFAC  370 

J=IC+:i  QFAC  380 

10  J=J-1  OFAC  390 

IF(J-1)420,420,20  CFAC  400 

20  IF(C( J)  )30, 10,30  QFAC  410 

C  QFAC  420 

C         NORMALIZATION  OF  REMAINING  COEFFICIENTS  OFAC  430 

30  A=C(J)  QFAC  440 

IF( A-l. 00)40, 60, 40  QFAC  450 

40  DO    50  I<UJ  QFAC  460 

C(l  )*C(  1  )/A  QFAC  470 

CALL  OVERFL(N)  OFAC  480 

lF(N-2>420,50, 50  OFAC  490 

50  CONTINUE  OFAC  500 

DNORMxA  OFAC  510 

C  QFAC  520 

C         TEST  ON  NECESSITY  OF  BAIRSTOW  ITERATION  QFAC  530 

60  I F ( J-3)430,380,70  QFAC  540 

C  QFAC  550 

C         PREPARE  BAIRSTOW  ITERATION  QFAC  560 

70  EPS-l.D-14  OFAC  570 

EPSU1.0-6  QFAC  580 

L«=0  OFAC  590 

LL=0  CFAC  600 

01=011)  QFAC  610 

02=012)  QFAC  620 

0Q1=0.D0  OFAC  630 

002=0.00  QFAC  640 
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AA=C( 1)  QFAC  650 

BB=C(2)  QFAC  660 

C8=DA6S(AA)  OFAC  670 

CA=DABS(BB)  QFAC  680 

IF(CB-CA)80,.90,100  OFAC  690 

80  CC=CB*CB  OFAC  700 

CB=C8/CA                                                .  OFAC  710 

CA=1.00  QFAC  720 

GO    TO       110  QFAC    730 

90    CC=CA  +  CA  QFAC    740 

CA=1.00  OFAC  750 

CB=1.D0  OFAC  760 

GO  TO   110  OFAC  770 

100  CC=CA+CA  QFAC  780 

CA=CA/CB  OFAC  790 

CB=1.00  QFAC  800 

110  CD=CC».1D0  OFAC  810 

C  QFAC  820 

C         START  BAIRSTOW  ITERATION  OFAC  830 

C         PREPARE  NESTED  MULTIPLICATION  QFAC  o'.O 

120  A=O.DO  QFAC  850 

B=.A  QFAC  660 

A1=A  QFAC  870 

B1=A  QFAC  880 

I=J  OFAC  890 

QQQUQ1                                                   .  QFAC  900 

QQQ2=02  OFAC  910 

DQ1=HH  OFAC  920 

DQ2=H  QFAC. 930 

C  QFAC  940 

C         START  NESTED  MULTIPLICATION  OFAC  950 

130  H=-Ql»B-Q2»AtC( I )  OFAC  960 

CALL  OVERFLIN)  QPAC  970 

IF(N-2)440, 140 i 140  OFAC  980 

140  B=.A  QFAC  990 

A=H  QFACLOOO 

1=1-1  QFAC1010 

IFU-1U80,  150,  160  0FAC1020 

150  H=0.D0  OFAC1030 

160  H=-Q1»B1-Q2»A1+H  CFAC1040 

CALL  OVERFLIN)  0FAC1050 

lF(N-2)440, 170, 170  QFAC1050 

170  C1=B1  CFAC1070 

61=A1  QFAC1C30 

A1  =  H  OFAC1090 

GO  TO   130  QFACUOO 

C         END  OF  NESTED  MULTIPLICATION  QFAC1110 

C  QFAC1120 

C         TEST  ON  SATISFACTORY  ACCURACY  0FAC1130 

180  H=CA»DABS(A)+CB»DABS(B)  CFAC1140 

IFILL) 190,190,390  OFAC1150 

190  L=L+:i  0FACU60 

IF(DABS( A)-EPS»DABS(CI 1) ) ) 200, 200, 210  QFAC 1170 

200  IF(DABS(3)-EPS»DABS(C(2)) 1390,390,210  QFAC1180 

C  QFAC 11 90 

C          TEST  ON  LINEAR  REMAINDER  OF  MINIMUM  NORM  QFACUOO 

210  IF(H-CC)220r220,230  QFAC1210 

220  AAxA  CFAC1220 

08=B  OFAC1230 

CC=H  0FAC1240 

QQ1=QI  0FAC1250 

QQ2=Q2  0FAC1260 

C  QFAC1270 

C                       TEST    ON    LAST     ITERATION    STEP  QFAC1280 

230    IFt L-LIK1280, 280,240  QFAC1290 

r  QFACUOO 
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C         TEST  ON  RESTART  OF  BAIRSTOW  ITERATION  WITH  ZERO  INITIAL  GUESS   OFAC1310 

240  1F(H-CD)450,.45D.250  CFAC1333 

250  1F(0(  1)  1270^260,270  QFAC1340 

260  IF(C(2> 1270,440,270  QFAC135D 

270  Olll^O.OO  .  QFAC1360 

0<2)=»O.DO  .  CFAC1370 

GO  TO    70  QFAC1380 

C  UFAC1390 

C         PERFORM  ITERATION  STEP      „....  QFAC1400 

280  HH'DHAXUDADSI  A  1  )  ,  DABS  (  B  1  >  ,  DAD  S  (  C  1  )  )  QFAC1410 

1F(HH)440,  440,290  .   •  CFAC1420 

290  A1  =  A1/HH  (JFAC1433 

B1  =  B1/HH  QFAC1440 

C1  =  C1/HH  QFAC1453 

H=A1»C1-B1»B1  QFAC1460 

IF(H)300,440,300  0FAC1473 

300  A=A/HH  QFAC1480 

B=B/HH  0FAC1493 

HH=( B*Al-A»81)/H  QFAC1500 

H=( A«C1-B»B1)/H  QFAC151D 

Q1=.Q1  +  HH  QFAC1523 

Q2=Q2+H  '  0FAC1530 

C         END  OF  ITERATION  STEP  QFAC1540 

C        TEST  ON  SATISFACTORY  RELATIVE  ERROR  OF  ITERATED  VALUES  %^c\Hl 

IF(DABS(HH)-EPS»DABS(Q1))310,.310,330  •  OFAC1570 

310  IF(DABS(H  -EPS.DABS(Q2))320,.320,330  GFAC15B0 

320  LL=1  0FAC159J 

GO  TO   120  .  QFAC1600 

C         TEST  ON  DECREASING  RELATIVE  ERRORS  QFici620 

330  IF(L-1)120, 120,340  OFAC163D 

340  IF(DABS(HH)-EPS1»DABSIQ1))350,350,120  oFAfl   0 

350  IFIDABS(H)-EPS1.DABS(Q2))360,.360,120  J**" 

360  IF(DABS(QQ01>HH)-DA3SIQI«DU1>.)370,460,460  CFAC1660 

370  IFtDABS(QQa2.H)-DABS(Q2.DQ2))120,460,460  OFAC167D 

C         END  OF  BAIRSTOW  ITERATION  (JFAC1680 

C         EXIT  IN  CASE  OF  QUADRATIC  POLYNOMIAL  2fAC1700 

380  0(1)=C(1)  0FAC1710 

0'2'cCt2»  0FAC1720 

Q(3>=0.00  QFAC1730 

Q(4)=0.D0  QFAC1740 

IC=J~2  0FAC175D 

RETURN  UFAC1760 

C         EXIT  IN  CASE  OF  SUFFICIENT  ACCURACY  SfAC1783 

390  Q|l)xQl  0FAC1793 

0<2  *02  CFACICOO 

Q,31cA  QFAC1810 


Q(4> 

DIVIDE  OUT  QUADRATIC  FACTOR 


CFAC1820 
QF AC  18 30 

B  -C   C-l)  CFAC1840 

.-  1  QFAC185D 


J=IC-2 
IC=lC-2 

400  C1=A1 

A1=B1-C1»Q2 

B1=C( JJ-C1»Q1 

C( J)=C1 

J=J-1 

IF(  J-D400,  410,400 
410  C(J)=Al 

RETURN 


CFAC1S50 
UFACie7D 
CFAC188D 
0FAC1890 
CFAC1900 
OF  AC  1910 
QFAC1920 
QFAC1930 
QFAC1940 
OFAC1950 


C         ERROR  EXIT  IN  CASE  OF  ZERO  OR  CONSTANT  POLYNOMIAL  QFAC196D 


144 


420  IER=-1  QFAC107J 

RETURN  CFAC1980 

C  QFAC1990 

C         ERROR  EXIT  IN  CASE  OF  LINEAR  POLYNOMIAL                          GFAC2000 

430  IER=-2  QFAC2010 

RETURN  QFAC2020 

C  0FAC2030 

C         ERROR  EXIT  IN  CASE  OF  NONREFINED  QUADRATIC  FACTOR               QFAC2040 

440  IER=-3  0FAC2050 

GO  TO   460  QFAC2060 

C  QFAC2070 

C         ERROR  EXIT  IN  CASE  OF  UNSATISFACTORY  ACCURACY                   QFAC2080 

450  IER=il  QFAC2090 

460  Q(1)«=Q«1  QFAC2100 

Q(2)=UQ2  UFAC2113 

IC=J-1  (JFAC2120 

Q(3)=AA  QFAC2130 

Q(4)=BB  QFAC2140 

RETURN  QFAC2150 

END  QFAC2i60 


SUBROUTINE  MAGSQ<  X  ,  Z..N  1  MGSQ   10 

C      SUBROUTINE  MAGSQ  COMPUTES  THE  MAGNITUDE  SQUARE  POLYNOMI AL ,  I  .  E .     MiiSQ   20 

C  .          MGSQ   30 

C                      2            2  •                         MGSQ   40 

C                    Z(W  )  =   X(JW)    =  X(JW)X(JW)  MGSQ   50 

C  MGSQ   60 

C      X  =  INPUT  POLYNOMIAL  MGSO   70 

C  2IN-1)           MGSQ   80 

C      Z  =  MAGNITUDE  SQUARE  POLYNOMIAL  Z(N)  =  COEF.  OF  W                 MGSQ   90' 

C      N  =  NUMBER  OF  COEFICISNTS  (ORDER  ♦  1)  MGSQ  100 

C  J.M.  ELDER   4/20/71     MGSQ  110 

OOUBLE  PRECISION  X,Z  .                    MGSQ  120 

DIMENSION  X(l),Zll)  MGSQ  130 

DO    10  1  =  1, N  MGSQ  140 

10  Zt I (=0.0  MGSQ  150 

DO    50  1  =  1, N  MGSQ  160 

DO    50  J  =  1,.N  MGSQ  170 

K=I+J  MGSQ  180 

C      TEST  FOR  ODD  K  (TERM  OF  1*2)... KEEP  ONLY  EVEN  K  MGSQ  190' 

FK=FLOAT(K) /2.0+1.0E-4  MGSQ  200 

IF(K-2»INT( FK) (50, 20,50  MGSQ  210 

20  K=INT(FK)  M^SUi  220 

C      DETERMINE  SIGN  OF  SUM  T ERM .. .CHANGE  SIGN  IF  (I+3J)/2  =  K+J  ODD     MGSO  230 

FK=FLOAT(K*J)/2.0+1.0E-4  MbSQ  240 

IF(2»II'IT(FK).-K-J)30,40,30  MGSQ  250 

30  ZU)=-X(  I  )»X(J)+Z(K)  MGSQ  260 

GO  TO    50  MGSQ  270 

40  Z(K)-Xm*X(J)+Z(K)  NliSQ  2d0 

50  CONTINUE  MGSQ  290 

RETURN  MGSQ  300 

END  MGSQ  310 
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DOUBLE  PRECISION  PHI ,PHIK , PS  I , F ACTR, VAL , Dl , D2 

DIMENSION  PHK  20),  PHIM20)  ,  PSI  (20) ,F ACTR (20) ,A2K(20) 

COMMON  A2K, PHI , PHIK.D1.D2 

DATA  I0LK/3H    /.IP/3HP   /.IPK/3HPK  /, I END/3HEND/ , I  GO/ 3H« 30/ 
10  FORMAT  I 1A3, 12,025.16) 

20  FORMAT (//11X.26HCHARACTERISTIC  POL YNOM I AL S / 39X , 1HN/2 2 X  , 
117HC0EFFICIENTS  OF  S/6X  ,  1HN  ,.5X  ,  14H0PE.N  LOOP  (  PH  I  )  ,  4  X  , 
217HCLOSED  LOOP(PHIK)) 
30  FORMAT  (  I7.1P2LU9.7  ) 

40  FORKAT(//5X»-10HDIAGO'JAL  Q/ 5X..1  7HDI  AGONAL  TERMS. . . / T22X , 1 PD1 9. 12 ) ) 
50  F0RMAT(//5X,.11HUNIT  RANK  Q/ 5X  ,.25HA  VER  AGE  ERRUR  IN  FACTOR  =, 
11P011.3/20X,6HFACT0R,20X,3HPSI  ,  21X  ,  5HCHECK  /  (  10  X  ,  1P3D2  5  .  12  )  ) 
60  FORMAT (5X,50H«««»»  POLYNOMIAL  HAS  HIGHEST  POWER  TERM  ZERO  •!»♦•) 
70  FORMAT!  5X,46H»»»««  INPUT  CONTAINS  UNI  DENTI F  1 6L6  LA3EL  •••••) 
80  F0RMAT(5X,47H»»»»»  SPECTRAL  FACTORIZATION  UNSUCCESSFUL  •••••) 
90  FORMAT(5X,35H»»"»  SYSTEM  PAIR  NON-OPTIMAL  •••*•) 
100  FORMATdH  ,A3,  I2.D25.16) 
110  FORMAT  KH  PHI) 
120  F0RMAU5H  PHIK) 

NDIM=4 

PHI(1)=5.0 

PHI(2)=9.0 

PHI I3)=5.0 

PHI(4)=1.0 

PHlK(l)-26.0 

PHI K( 2) =25.0 

PH1KI3)=8.0 

PHI K( 4  1  =  1.0 

GO  TO   280 
130  DO   1-0  1=1,20 

PHI (  I  , =0.000 
143  PHIKt I  1=0.000 
150  JMP=-1 

NDIM*0 

OECODE  INPUT 
160  READI5,   10 ) ICHAft, I, VAL 

HRITE(6,.  100)1CHAR,  I..VAL 

1F(  ICHAR-IBLK.1170,210,  170 
170  IF( ICHAR-IP)1S0,240,  180 
180  IF(  ICHAR-IPK1190, 220,-190 
190  IF( ICHAR-IG0)200, 280,200 
200  1F( ICHAR-IEN01420, 430,420 
210  IF(JMP) 160,250,230 
220  JMP=1 
230  PHIKI I )=VAL 

WRITEI6,.  120) 

GO  TO   260 
240  JMP  =  0 
250  PHI(I)=VAL 

WRITE16,.  110) 
260  IFINOIM-I )270,  160,  160 
270  NDIM=I 

GO  TO   160 

NORMALIZE  POLYNOMIALS  SO  THAT  HIGHEST  POWER  TERM  IS  ONE 
280  N0RD=N0IM-1 

01=PHI INDIHl 

D2=PHIK(NDIM) 

IF(DABS(D1 1-1.00-40)410,410, 290 
2  90  IF (DA0SID2 1-1.00-40)4 10, 410, 300 
300  00   310  I=1,N0RD 

PHI  ( 1  )=PHI (  I  1/01 
310  PHIKt  I  )=.PHIK(  I  1/02 

PHI (NDIM)=1 .ODO 

PHIK(NOIM)= 1.000 


OPTI 

10 

OPTI 

20 

UPTI 

30 

OPTI 

40 

OPTI 

50 

OPTI 

60 

OPTI 

70 

OPTI 

80 

OPTI 

93 

CPTI 

100 

OPTI 

110 

OPTI 

120 

OPTI 

130 

OPTI 

140 

OPTI 

150 

OPTI 

160 

OPTI 

170 

OPTI 

180 

OPTI 

190 

OPTI 

200 

OPTI 

210 

OPTI 

220 

OPTI 

230 

OPTI 

240 

OPTI 

250 

OPTI 

260 

OPTI 

270 

OPTI 

280 

OPTI 

290 

OPTI 

300 

OPTI 

313 

OPTI 

320 

OPTI 

333 

OPTI 

340 

OPTI 

350 

OPTI 

360 

OPTI 

370 

OPTI 

380 

OPTI 

390 

OPTI 

400 

OPTI 

413 

OPTI 

420 

OPTI 

433 

OPTI 

440 

OPTI 

453 

OPTI 

463 

OPTI 

470 

OPTI 

480 

OPTI 

493 

OPTI 

500 

OPTI 

513 

OPTI 

520 

OPTI 

533 

OPTI 

543 

OPTI 

550 

OPTI 

563 

OPTI 

573 

OPTI 

580 

UPTI 

593 

OPTI 

600 

OPTI 

613 

OPTI 

620 

OPTI 

630 

OPTI 

640 
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320 


330 


KRITE16,.   20) 

00   320  I  =  1,.NDIM 

J=I-1 

WRITEI6,   30)J,PHI(  Il.PHIKU) 

COMPUTE  MAGNITUDE  SCUA3E  POLYNOMIALS 
CALL  MAGSUI PHI .PSI.NDIM) 
CALL  MAGSQt PHIK, PHI, NDIM) 


3*0 


360 
370 
380 

390 
*C0 
410 
420 
430 


PSI=  PHI 


PHI  =  PHIK 


COMPUTE  DIFFERENCES  OF  MAGNITUDE  SQUARE  POLYNOMIALS 
00   340  I  =  WNDIM 
PSII  I)  =  PHI(  I  J-PSK  I) 
A2KII)=SNGL(PSI (  I)  ) 

2 
IS  THE  CLOSED-LOOP  SYSTEM  OPTIMAL, I.E.  IS  PSItW  )  POSITIVE  SEMIO 
ITEST=NCNNEGIA2K,N0RD, 1) 
IF( ITEST-5)39D,350i360 
DIAGONAL  0 

WRITE16,.   40)1  PSK  I  1,1  =  1, NORD) 
RANK  ONE  Q 

CALL  SPFCTR(PSI,FACTR,NORD, IER) 
IF(  IER)400,370,400 
CHECK  FACTORIZATION 
CALL  MAGSOt  FACTR.PHUNORD) 
01=0.000 

DO   380  1  =  1, .NORD 
D1=D1+DABS<  PSH  I  )-PHI(  I  )  ) 
D1=01/D9LE( FLOAT (NORD) ) 

WRITE! 6,   50)D1, ( FACTR(I) ,PSI ( I), PHI (I  1,1  =  1, NORD) 
GO  TO   130 
WRITEI6,.   90) 
GO  TO   130 
WRITEC6*   80) 
GO  TO   130 
HRITE16,   60) 
GO  TO   130 
WRITEI6,   70) 
GO  TO   130 
STOP 
END 


OPTI  650 

OPTI  660 

OPTI  670 

OPTI  660 

GPTI  693 


3PTI 


'00 


GPTI  710 

CPTI  720 

CPTI  730 

OPTI  740 

OPTI  750 

CPTI  760 

OPTI  7  70 

OPTI  780 

JPTI  793 

OPTI  SCO 

OPTI  610 

OPTI  620 

OPTI  630 

OPTI  BO 

OPTI  850 

OPTI  660 

OPTI  870 

OPTI  830 

OPTI  890 

OPTI  900 

OPTI  910 

OPTI  920 

OPTI  930 

OPTI  9*0 

UPTI  950 

OPTI  960 

CPTI  970 

OPTI  9cO 

OPTI  990 
OPTI 100  0 
OPTI1010 
OPTI 102? 
OPTI 1030 
OPTI 1040 


pnorpAi'  opt  i 


CHARACTERISTIC    PCLYI.'Cf.l  ALS 
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H 


COEFFICIENTS  OF  S 


OPFIJ    LOOPCPIM  ) 

2.U9GC000D  Olt 

*.I|128C00D  01. 

3.230U00PO  Oti 

1.279COO0n  OU 

2.98200000  05 

fc.ioooooon  02 

3.1000000O  01 

1.00000000  00 


CLOSER    LOOPCPH   K) 

2.U9EOO0OD  Oil 

Ii.li2217870  OU 

3.2UE99?EP  OU 

1.291G180D  OU 

3.027E270D  03 

U.1937EU5P  02 

3.1E709U1D  01 

l.OOOOOOOD  00 


I  R0lf#  I 


NORMALIZED  ROUTt:  ARRAY 
NORMALIZED  ROI.'S 


•  R  1 

1 

t 

U.29E-01 

l.OOE  00 

5.33E-01 

-1.U3E-01 

-1 

C7E-01 

0 

0 

1.19E-C2 

2 

♦ 

S.1UE-01 

l.OOE  00 

U.67F.-01 

-8.57E-02 

-C 

E7E-12 

0 

0 

0.0 

3 

♦ 

8.57E-01 

l.OOE  00 

-3.G7E-01 

-5.71E-01 

0 

0 

6 

12E-02 

t> 

♦ 

5.82E-01 

l.OOE  00 

3.7UE-01 

-9.70E-02 

-5 

35F.-02 

0 

0 

5 

- 

-5.1UE-01 

-l.OOE  00 

-U.E7E-01 

8.57E-02 

6 

E7E-02 

6 

- 

-8.57E-01 

-l.OOE  00 

-1.CUE-1G 

1.U3E-01 

0 

0 

•  R  2 

7 

- 

-8.57E-01 

-l.OOE  00 

-2.U3E-16 

1.U3E-01 

8 

- 

-l.OOE  00 

-7.7EE-01 

-9.UUE-17 

0.0 

9 

- 

-l.OOE  00 

-U.8f.E-16 

U.29E-01 

10 

- 

-l.OOE  00 

-5.51E-01 

0.0 

11 

+ 

l.OOF.  00 

7.78E-01 

12 

+ 

l.OOE  00 

0.0 

•R  3 

13 

* 

l.OOE  00 

1  ROUS  PRECEEPINP  ZERO  ROI.'S 
(ZERO  P.r»i:S  HAVE  BEEII  RESOLVED  EY  D-l  FFERENT I  ATI  IIG  PRECEEDI  l.'C  POl') 


PSKl.1  )  IS POSITIVE  SEP!  DEFINITE,  POLYNOMIAL  PASSES  SUFFICIENCY  TEST 


UNIT  RANK  Q 
AVERAGE  ERROR 


IN   FACTOR    '     6.15 

FACTOR 
0000000000000    00 
UGU0161S13ED    00 

oooooooonooco  oo 

039230'iE!iSUin  CI 

loooooooooccn  01 

92J.20323027CO  00 
OOOOOOOOOOOCP  00 


PS  I 

.OOOOOOOCOOOOP  00 
0 

LOCOOCOOOOOPP  01 
200000000COCP  01 
9000POOOOOOOP  01 

dooocorooooo?  01 

6000000000000  01 


CHFCK 
l.OOCOOOCCOCOCP  00 
li.UU089?09F5010-lS 
-l.UOPOOOroCCOP^  CI 

1. 2rooooroccoon  01 

U.9PC0n0P0P00CP  01 

-8.upppococrcr-cp  ci 

3.6OC000000000P  01 
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